Algoritmo para la estimación de la turbidez sobre el Río Paraná

Autor/a
Afiliación

Víctor Gustavo Gómez

GISTAQ (UTN-FRRe)

Fecha de publicación

13 de octubre de 2025

Resumen

Este sitio web contiene información sobre la estimación de la turbidez por teledetección en la cuenca media del Río Paraná. La turbidez es uno parámetros de interés dentro proyecto Estimar indicadores de calidad de agua en la cuenca media del río Paraná para el desarrollo de un algoritmo mediante técnicas de teledetección satelital (MSECRE0008604), desarrollado por el Grupo de Investigación Sobre Temas Ambientales y Químicos (GISTAQ) de la Universidad Tecnológica Nacional Facultad Regional Resistencia (UTN-FRRe).

Se utilizarán imágenes del satélite Sentinel-2 con corrección automática, de las cuales se obtiene la reflectancia de superficie del agua. Se buscará la relación entre la reflectancia y la turbidez por métodos de regresión tradionales y machine learning. Una vez obtenido el algoritmo que relacione ambas propiedades, se desarrollaran mapas de distribución espacial.

Palabras clave

GISTAQ, UTN, FRRe, Algoritmo, Turbidez, Machine learning, Teledetección

1 Turbidez

La turbidez se refiere a la opacidad o falta de claridad en un líquido provocada por la presencia de partículas suspendidas. Este fenómeno es un indicador clave en el monitoreo de la calidad del agua y su influencia en diferentes ecosistemas es significativa.

La turbidez es un indicador de la calidad del agua, reflejando la presencia de partículas en suspensión. Su medición es crucial para garantizar la potabilidad del agua y la salud de los ecosistemas acuáticos. Este fenómeno puede ser resultado de diversas causas, como la erosión del suelo, la actividad biológica y la contaminación. La comprensión de la turbidez y su impacto es esencial para la gestión de recursos hídricos y la protección del medio ambiente.

La turbidez viene determinada por la dispersión de la luz causada por la materia suspendida en el agua, se obtiene normalmente mediante un turbidímetro, que proporciona medidas en Nephelometric Turbidity Unit (NTU) y mide la dispersión de un rayo de luz en el agua a 90º de la luz incidente [1].

Muchas propiedades, como la clorofila-a (Chl-a), sólidos suspendidos totales (SST) y la materia orgánica disuelta coloreada (CDOM), se utilizan a menudo como indicadores del estado del agua. Estos constituyentes del agua a su vez son responsables de la turbidez.

Existe una fuerte correlación entre turbidez y sólidos suspendidos totales, por lo que se puede estimar SST a partir de la turbidez. Por lo general, es una relación directa, a mayor concentración de SST mayor turbidez.

Existe una relación inversa entre la Turbidez y la profundidad del disco de Secchi (a valores bajos de secchi mayor turbidez), por lo que también se puede estimar turbidez a partir de mediciones de disco de secchi.

1.1 Métodos tradicionales

Tabla 1: Características principales de algoritmos tradicionales para la estimación de turbidez.
Ecuación Bandas (nm) Métricas Aguas Plataforma Referencia
\(1.559e^{35.533 \cdot B03} \\ 1.879e^{37.745(B03 \cdot B5)/(B04+B12)}\) B03, B04, B05, B12 \(R^{2}\), RMSE, MAE Lago1 Sentinel-2 [2]
\(2677.2 \cdot B04^{1.856}\) B04 \(R^{2}\), RMSE, Bias Interiores variadas2 Landsat-8 [3]
\(969-1.5468 \cdot R_{1200nm}+2.07 \frac{B8A}{B02}\) B02, B8A, 1200nm IOA, SI, RMSE, MAE Río3 Landsat-8 [4]
\(y=-1.1+5.8 \frac{B02}{B04} \\ y=3.896-4.186 \frac{B02}{B03}\) B02, B03, B04 \(R^{2}\), RMSE Río4 Landsat-8 [5]
\(y=37661 \cdot B8A^{2}+1845 \cdot B8A <br> y=531.5- \frac{B04}{0.88}\) B04, B8A \(R^{2}\), RMSE, MAPE Estuario5 Pléiades [6]

Múltiples modelos (lineal, logaritmos, inversa, cuadrática, exponencial, potencial) y plataformas (Sentinel-2, Landsat-5 y Landsat-8) emplean el cociente de bandas B04/B03 [7].

Modelos de estimación a partir de Sentinel-2 y Landsat-8 utilizan regresiones lineales, cuadráticas y logarítmicas empleando B02, B03, B04, B01 (con menos apariciones) y cocientes entre éstas [8].

1.2 Métodos de aprendizaje automático

El aprendizaje automático es un subconjunto de la inteligencia artificial que permite que un sistema aprenda y mejore de forma autónoma, sin necesidad de una programación explícita, a través del análisis de grandes cantidades de datos. El aprendizaje automático permite que los sistemas informáticos se ajusten y mejoren continuamente a medida que acumulan más “experiencias”. Por lo tanto, el rendimiento de estos sistemas puede mejorar si se proporcionan conjuntos de datos más grandes y variados para su procesamiento.

Cuando se entrenan modelos de machine learning, cada conjunto de datos y cada modelo necesitan un conjunto diferente de “hiperparámetros”. Los hiperparámetros son variables de configuración externa que se utilizan para administrar el entrenamiento de modelos de machine learning. Controlan de forma directa la estructura, funciones y rendimiento de los modelos. Los hiperparámetros son los parámetros de un modelo de aprendizaje automático, que no se aprenden durante el entrenamiento, sino que se establecen antes de que comience.

El “ajuste de hiperparámetros” permite modificar el rendimiento del modelo para lograr resultados óptimos. Este proceso es una parte fundamental del machine learning. El ajuste de hiperparámetros puede ser manual o automático. A pesar de que el ajuste manual es lento y tedioso, permite entender mejor cómo afectan al modelo las ponderaciones de los hiperparámetros. El proceso de ajuste de hiperparámetros es iterativo, y debe probar diferentes combinaciones de parámetros y valores.

En el aprendizaje automático es importante utilizar técnicas de “validación cruzada” , de modo que el modelo no se centre únicamente en una única porción de sus datos. La validación cruzada o cross-validation es una técnica utilizada para evaluar los resultados de un análisis estadístico y garantizar que son independientes de la partición entre datos de entrenamiento y prueba. La idea básica de la validación cruzada es dividir los datos en conjuntos de entrenamiento y validación, y luego entrenar el modelo en el conjunto de entrenamiento y evaluar su rendimiento en el conjunto de validación. Este proceso se repite varias veces, con diferentes subconjuntos de los datos utilizados para el entrenamiento y la validación, y se calcula el rendimiento promedio.

En los procesos de machine learning supervisado se utilizan diversos algoritmos y técnicas de cálculo, generalmente calculados mediante el uso de programas como R o Python.

Dependiendo del tipo de datos que se usen para el entrenamiento, será de modelo de aprendizaje automático que se use. A grandes rasgos, existen tres tipos de modelos que se usan en el aprendizaje automático: aprendizaje supervisado , no supervisado y por refuerzo.

Consultando el trabajo de otros investigadores, se observa que utilizan principalmente el aprendizaje automático supervisado. Este tipo aprendizaje supervisado utiliza un conjunto de entrenamiento para enseñar a los modelos a producir el resultado deseado. Este conjunto de datos de entrenamiento incluye entradas y salidas correctas, que permiten al modelo aprender con el tiempo. El algoritmo mide su precisión a través de la función de pérdida, ajustando hasta que el error se haya minimizado lo suficiente.

Yang Zhe y otros, utilizaron como datos de entrada la reflectancia de superficie y datos de salida la turbidez, utilizaron los modelos SVR (support vector regression), random forest (RF) y eXtreme Gradiente Boostring (XGBoost). Los hiperparámetros de cada modelo se determinaron mediante una búsqueda en cuadrícula de validación cruzada en Scikit-Learn de Python [9].

Ma Yue y otros, utilizaron varios modelos de aprendizaje automático, usaron Python 3.7 tanto para la predicción de la turbidez del agua y la optimización de la los hiperparámetros [2].

Zhao y otros probaron 14 modelos de machine learning en un estanque de peces con un dispositivo de construction propia, de los cuales ETR, Bagging, RFR, and ABR son los que presentaron un mejor desempeño en la estimación de la turbidez. Los algoritmos se implementaron utilizando Python 3.6 y bibliotecas de aprendizaje scikit [10].

Tabla 2: Características principales de algoritmos de aprendizaje automático para la estimación de turbidez.
Modelo de machine learning Cuerpo de agua Métricas Plataforma Referencia
SVR, ELM ,BP ,CART ,GBT ,RF ,KNN Lagos RMSE, \(R^{2}\), MAE Sentinel-MSI [2]
eXtreme Gradient Boosting (XGBoost), support vector regression (SVR), random forest (RF) Lago RMSE, \(R^{2}\), MAPE Sentinel-2A/B y Landsat-8/9 [9]
linear regression (LR), ridge regression (RR), least absolute shrinkage and selection operator regression(LASSO), elastic net regression (ENR), k-nearest neighbor regression (KNN), Gaussian process regression (GPR), decision tree regression (DTR), support vector regression (SVR), multilayer perceptron regression (MLP), adaptive boosting regression (ABR), gradient boosting regression (GBR), bootstrap aggregating regression (Bagging), random forest regression (RFR), and extreme tree regression (ETR) Estanque de peces MAE, MRSE, MAPE, \(R^{2}\), RE, Acc Propia [10]

2 Procesamiento de datos

Para el procesamiento de los datos se utilizará la librería pandas de Python.

En el proyecto tenemos dos archivos .csv que contienen los datos:

-base_de_datos_lab.csv → contiene resultados de laboratorio

-base_de_datos_gis.csv → contiene datos espectrales

Importamos la librería pandas para usarla, la nombramos como “pd” para simplificar

import pandas as pd

Leemos los archivos .csv por separado y definimos dos DataFrame

Un DataFrame es basicamente una tabla, donde la información se organiza en filas y columnas. Los datos de la misma columna contienen el mismo tipo de datos, pandas agrega por defecto un “índice” que nos ayuda a identificar una fila en particular.

Con la función pd.read_csv le idicamos a pandas que queremos leer archivos .csv.

df1_lab = pd.read_csv(r"D:\GIT\Turbidez\datos\base_de_datos_lab.csv")
df2_gis = pd.read_csv(r"D:\GIT\Turbidez\datos\base_de_datos_gis_acolite.csv")

df1_lab DataFrame de datos provenientes del laboratorio.

df2_gis DataFrame de datos espectrales provenientes del sensor MSI de Sentinel-2.

Nota: Se debió colocar la “r” delante de la dirección para que lea los archivos.

Video de YouTube ¿Qué es un DataFrame?

Verificamos que los datos se han leído y se crearon correctamente ambos DataFrame por separado, con print

df1_lab.head()

df2_gis.head()
punto banda reflect fecha longitud latitud
0 1 B01 0.041367 2023-05-11 -58.868047 -27.464687
1 1 B02 0.046669 2023-05-11 -58.868047 -27.464687
2 1 B03 0.085630 2023-05-11 -58.868047 -27.464687
3 1 B04 0.128799 2023-05-11 -58.868047 -27.464687
4 1 B05 0.128704 2023-05-11 -58.868047 -27.464687

Nota: La primer columna (donde se ven los valores 0,1,2,3,4, …) es el índice que agrega pandas por defecto al DataFrame, esa columna no forma parte del csv original.

Combinamos ambos DataFrame para tener los datos en una única base de datos.

Lo defimos como df_combinado, para esto utilizamos la función pd.merge para realizar la combinación.

df_combinado = pd.merge(df1_lab, df2_gis, on=['latitud', 'longitud','fecha'], how='inner')

on=[‘latitud’, ‘longitud’,‘fecha’] Especifica las columnas por las cuales se unirán los dos DataFrames. En este caso, por coincidencias exactas en latitud, longitud y fecha.

how=‘inner’ Es el tipo de unión, significa que solo se conservarán las filas que tengan coincidencias en ambos DataFrames en las columnas mencionadas.

Verificamos que la combinacion se haya realizado correctamente

df_combinado.head()
fecha longitud latitud param valor punto banda reflect
0 2023-05-11 -58.868047 -27.464687 ph 6.8 1 B01 0.041367
1 2023-05-11 -58.868047 -27.464687 ph 6.8 1 B02 0.046669
2 2023-05-11 -58.868047 -27.464687 ph 6.8 1 B03 0.085630
3 2023-05-11 -58.868047 -27.464687 ph 6.8 1 B04 0.128799
4 2023-05-11 -58.868047 -27.464687 ph 6.8 1 B05 0.128704

Filtramos la turbidez del DataFrame

Del DataFrame combinado, solo nos interesa las filas que contengan los valores de turbidez, ya que es la propiedad de estudio en este sitio web. Por ello nos quedamos con las filas en donde la columna param sea igual a turb.

Creamos un nuevo DataFrame a partir de df_combinado.

df_turbidez = df_combinado[(df_combinado['param'] == 'turb')]

Verificamos

df_turbidez.head()
fecha longitud latitud param valor punto banda reflect
33 2023-05-11 -58.868047 -27.464687 turb 185.5 1 B01 0.041367
34 2023-05-11 -58.868047 -27.464687 turb 185.5 1 B02 0.046669
35 2023-05-11 -58.868047 -27.464687 turb 185.5 1 B03 0.085630
36 2023-05-11 -58.868047 -27.464687 turb 185.5 1 B04 0.128799
37 2023-05-11 -58.868047 -27.464687 turb 185.5 1 B05 0.128704

Nota: Se han eliminado las filas de ph, cond, sol_sus, etc. Se conserva el índice del DataFrame original.

Eliminamos las columnas que no nos interesan

Creamos un nuevo DataFrame a partir de df_turbidez. Con la función .drop especificamos que columnas queremos eliminar.

df_turbidez_banda = df_turbidez.drop(columns=['longitud','latitud','punto','fecha','param'])

df_turbidez_banda.head()
valor banda reflect
33 185.5 B01 0.041367
34 185.5 B02 0.046669
35 185.5 B03 0.085630
36 185.5 B04 0.128799
37 185.5 B05 0.128704

Nota: Solo nos quedamos con las columnas valor,banda y reflect.

Cambiamos el nombre de la columna “valor” por el de “turbidez” Como los valores de esa columna son de turbidez, directamente cambiamos el nombre de la columna con la función .rename

df_turbidez_banda.rename(columns={'valor': 'turbidez'}, inplace=True)

df_turbidez_banda.head()
turbidez banda reflect
33 185.5 B01 0.041367
34 185.5 B02 0.046669
35 185.5 B03 0.085630
36 185.5 B04 0.128799
37 185.5 B05 0.128704

Creamos la tabla final

Usamos la función .pivot_table para que las columnas sean la turbidez y las distintas bandas de Sentinel-2 (B01, B02 , B03…)

df_final = df_turbidez_banda.pivot_table(
    index='turbidez',
    columns='banda',
    values='reflect',  
    )

¿Qué significa cada término en el argumento de la función pivot_table?

index=‘turbidez’ → Las filas serán los valores únicos de ‘turbidez’

columns=‘banda’ → Las columnas serán los valores únicos de ‘banda’ (B01, B02…)

values=‘reflect’ → El contenido de la tabla será lo que haya en la columna ‘reflect’

Creamos archivo final con los valores de turbidez y reflectancia de cada banda

-Forma 1 Creamos un .csv con los datos que nos interesan con la función .to_csv

 df_final.to_csv('Turbidez_FINAL.csv', index=True)

IMPORTANTE: index=True , debe ser así para que el índice que definimos al pivotar sea una columna visible en el archivo csv.

Recordemos que la columna index es solo visible solo por la librería pandas, Si guardamos con index=False, se omite y no se guarda en el csv.

-Forma 2 Si decidimos poner index=False tenemos que usar una función adicional antes de exportar, debido a que la turbidez está en el índice, no como una columna.

Por lo tanto, luego de hacer el pivot se debe agregar una línea de código.

df_final = df_turbidez_banda.pivot_table(
    index='turbidez',
    columns='banda',
    values='reflect',  
    )

df_final = df_final.reset_index() #Línea de código que se debe agregar

reset_index() convierte el índice definido previamente como “turbidez” en una columna normal y reemplaza el índice por uno numérico estándar (0, 1, 2…). que es el predeterminado por pandas.

Finalmente exportamos como archivo csv.

df_final.to_csv('Turbidez_FINAL2.csv', index=False)

Verificación tabla final

import pandas as pd

df = pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')

df.head()
turbidez B01 B02 B03 B04 B05 B06 B07 B08 B11 B12 B8A
0 2.0 0.030264 0.035694 0.052761 0.038619 0.028929 0.017878 0.018247 0.014811 0.008830 0.006162 0.015677
1 3.0 0.018744 0.027450 0.044980 0.028040 0.018781 0.007693 0.007571 0.007623 0.001773 0.000112 0.003305
2 4.0 0.029141 0.036459 0.054631 0.039793 0.029767 0.022395 0.022966 0.019097 0.009430 0.007976 0.020290
3 4.5 0.022306 0.027934 0.045472 0.028727 0.020835 0.014462 0.017479 0.013988 0.000088 0.000361 0.013232
4 5.0 0.029902 0.030043 0.047235 0.033111 0.024060 0.011116 0.012519 0.007006 0.001982 0.000181 0.009074

3 Prueba de modelo lineal

Regresión lineal con método de mínimos cuadrados

Reemplamos el código del tutorial de sklearn por nuestros datos.

Debemos importar las funciones necesarias para:

1- Leer el archivo .csv creado previamente , que contiene los datos de turbidez y los valores de reflectancia para cada banda.

2- Para usar el método de mínimos cuadrados

Importamos las funciones

import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt

pandas → para leer los datos

train_test_split → para dividir los datos en entrenamiento y validación

LinearRegression → para crear el modelo de regresión lineal

mean_squared_error y r2_score → Para medir el desempeño del modelo (RMSE y R2)

matplotlib.pyplot → para realizar gráficos

Leemos los datos de interés y los dividimos en entrenamiento y validación.

Lectura de datos, lo hacemos con pd.read_csv

df = pd.read_csv(r"D:\GIT\Turbidez\Turbidez_FINAL.csv")

X = df[['B05']]  
y = df['turbidez']

Definimos:

df, un DataFrame que contiene los datos de turbidez y reflectancia;

X para que sea los valores de reflectancia de una banda en particular;

y para que sea la turbidez.

Dividimos en datos para el entrenamiento y validación

X_train, X_test, y_train, y_test = train_test_split(
  X, y, test_size=0.25, shuffle=True, random_state=42 
)

Conjunto de datos para el entrenamiento X_train y X_test

Conjunto de datos para la validación y_train y y_test

Creamos el modelo de regresión

regressor = LinearRegression().fit(X_train, y_train)

LinearRegression() Crea el modelo de regresión lineal llamado “regressor”. Este modelo, se entrena con la función .fit a partir de los datos de entrenamiento (X_train, y_train)

Evaluamos el modelo generado a partir de las métricas de desempeño.

y_pred = regressor.predict(X_test)
p_rmse = mean_squared_error(y_test, y_pred)
p_r2 = r2_score(y_test, y_pred)

Con la función print visualizamos los valores de las métricas de desempeño

print("RMSE", p_rmse) 
print("R2", p_r2)
RMSE 1151.59791674541
R2 0.27485564458517464

Para ver la ecuación del modelo de regresión

# Coeficientes (pendientes)
coef = regressor.coef_

# Intercepto (ordenada al origen)
intercept = regressor.intercept_

print(f"La ecuación es: y = {coef[0]:.3f}x + {intercept:.3f}")
La ecuación es: y = 1316.107x + -41.516

Visualizamos los resultados comparando el conjunto de entrenamiento y validación.

Código 
fig, ax = plt.subplots(ncols=2, figsize=(10, 5), sharex=True, sharey=True)

#Gráfico de entrenamiento
ax[0].plot(
    X_train,
    regressor.predict(X_train),
    linewidth=3,
    color="#17A77E",
    label="Modelo",
)

ax[0].scatter(X_train, y_train, label="Entrenamiento", color = "#9D50A6", alpha = .6)
ax[0].set(xlabel="Banda 5", ylabel="Turbidez", title="Conjunto de entrenamiento")
ax[0].legend()

#Gráfico de validación
ax[1].plot(X_test, y_pred, linewidth=3, color="#17A77E", label="Modelo")
ax[1].scatter(X_test, y_test, label="Validación", color = "#9D50A6", alpha = .6)
ax[1].set(xlabel="Banda 5", ylabel="Turbidez", title="Conjunto de validación")
ax[1].legend()

#Ecuación de la recta
coef = regressor.coef_[0]
intercept = regressor.intercept_
equation = f"y = {coef:.2f}x + {intercept:.2f}"
# Mostrar la ecuación en ambos subgráficos (opcionalmente, puedes usar solo uno)
for a in ax:
    a.text(0.05, 0.95, equation, transform=a.transAxes,
           fontsize=10, verticalalignment='top',
           bbox=dict(boxstyle="round", facecolor="white", alpha=0.7))

fig.suptitle("Regresión lineal")

plt.show()

4 Pruebas de correlación

Para ser más rigurosos, agregamos más etapas durante el entrenamiento de nuestro modelo lineal.

En esta etapa nos interesa estimar la turbidez por medio de modelos lineales sencillos, para encontrar mejores modelos podemos aplicar transformaciones logarítmicas con el objetivo principal de linealizar relaciones no lineales.

Primeramente se probará un modelo sin aplicar ninguna transformaciones en las variables.

En las siguientes secciones, se aplicará logaritmo en la turbidez, en las bandas espectrales y en en ambas, con el fin de evaluar en rendimiento de los modelos y seleccionar los mejores, según sus métricas de desempeño.

Se utilizará el coeficiente de correlación lineal r, para evaluar la relación lineal entre la turdidez y las distintas bandas (con y sin transformación logarítmica).

Esta medida indica cuanto se relacionan dos variables, puede tomar valores desde -1 a +1:

• +1 correlación lineal perfecta positiva.

• 0 sin correlación.

• -1 correlación lineal perfecta negativa.

4.1 Prueba de correlación C1: (sin aplicar logaritmo)

Importamos pandas para leer los datos

Código 
import pandas as pd

Datos= pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')
#print (Datos.head())

Calculamos coeficiente de correlacion lineal “r” entre la turbidez y cada banda con la función .corr de pandas.

Código 
bandas = [col for col in Datos.columns if col.startswith('B')]

correlaciones = {}

for banda in bandas:
    r = Datos['turbidez'].corr(Datos[banda])
    correlaciones[banda] = r

#Creamos un Data Frame.
df_correlaciones1 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones1 = df_correlaciones1.sort_values(by='r', ascending=False).reset_index(drop=True)

df_correlaciones1['Banda'] = "turb vs " + df_correlaciones1['Banda'].astype(str)
df_correlaciones1.rename(columns={'Banda': 'Combinación 1'}, inplace=True)

df_correlaciones1.to_csv(r'D:\GIT\Turbidez\Datos creados\Correlaciones\C1_turb_vs_banda.csv', index=False)


#Creo archivo csv para usarlo para entrenar modelos
Datos.to_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv", index=False)

df_correlaciones1
Combinación 1 r
0 turb vs B05 0.887299
1 turb vs B06 0.876827
2 turb vs B08 0.869149
3 turb vs B07 0.867741
4 turb vs B04 0.827019
5 turb vs B8A 0.788088
6 turb vs B03 0.673358
7 turb vs B02 0.614315
8 turb vs B01 0.611305
9 turb vs B12 0.429200
10 turb vs B11 0.409447

Gráfica Turbidez en funcion de las distintas bandas

Código 
import pandas as pd
import matplotlib.pyplot as plt

# Cargar los datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")

# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'turbidez']

# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'turbidez']

# Cantidad de bandas
n = len(bandas)

# Calcular las filas y columnas para el grid de subplots
filas = (n + 2) // 3  # hasta 3 columnas por fila
columnas = 3

# Crear figura y ejes
fig, axes = plt.subplots(filas, columnas, figsize=(12, filas * 3))
axes = axes.flatten()

# Graficar cada banda
for i, banda in enumerate(bandas):
    ax = axes[i]
    ax.scatter(df[banda], df['turbidez'], alpha=0.6, s=10)
    ax.set_title(f'Turbidez vs {banda}', fontsize=10)
    ax.set_xlabel(banda)
    ax.set_ylabel('Turbidez')
    ax.grid(True)

# Ocultar los subplots vacíos si hay
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

# Ajustar espacio
plt.tight_layout()
plt.show()

Las mejores bandas para predecir turbidez (por correlación lineal) son:

B05 > B06 > B08 > B07 > B04

4.2 Prueba de correlación C2: logaritmo en la turbidez

Importamos nunpy para operar con funciones matemáticas

Creamos un nuevo DataFrame para aplicarle el logaritmo a la colomna turbidez.

Código 
import numpy as np

Datos_turb_log = pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')
Datos_turb_log['turbidez'] = np.log(Datos_turb_log['turbidez'])

#Cambio el nombre la columna "turbidez" luego de aplicar el logaritmo
Datos_turb_log = Datos_turb_log.rename(columns={'turbidez': 'ln_turbidez'})

#Creo archivo csv para usarlo para entrenar modelos
Datos_turb_log.to_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C2_ln_turb_banda.csv", index=False)
#print(Datos_turb_log.head())

Calculamos r entre el ln_turb y cada banda, con la función .corr de pandas.

Código 
bandas = [col for col in Datos_turb_log.columns if col.startswith('B')]

correlaciones = {}

for banda in bandas:
    r = Datos_turb_log['ln_turbidez'].corr(Datos_turb_log[banda])
    correlaciones[banda] = r

#Creamos un Data Frame.
df_correlaciones2 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones2 = df_correlaciones2.sort_values(by='r', ascending=False).reset_index(drop=True)

df_correlaciones2['Banda'] = "ln_turb vs " + df_correlaciones2['Banda'].astype(str)
df_correlaciones2.rename(columns={'Banda': 'Combinación 2'}, inplace=True)

df_correlaciones2.to_csv(r'D:\GIT\Turbidez\Datos creados\Correlaciones\C2_ln_turb_vs_banda.csv', index=False)

df_correlaciones2
Combinación 2 r
0 ln_turb vs B05 0.813367
1 ln_turb vs B04 0.752422
2 ln_turb vs B06 0.714446
3 ln_turb vs B07 0.707941
4 ln_turb vs B08 0.703117
5 ln_turb vs B8A 0.603662
6 ln_turb vs B03 0.514314
7 ln_turb vs B02 0.429268
8 ln_turb vs B01 0.391705
9 ln_turb vs B12 0.242050
10 ln_turb vs B11 0.222510

Gráfica

Código 
import pandas as pd
import matplotlib.pyplot as plt

# Cargar los datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C2_ln_turb_banda.csv")

# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'ln_turbidez']

# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'ln_turbidez']

# Cantidad de bandas
n = len(bandas)

# Calcular las filas y columnas para el grid de subplots
filas = (n + 2) // 3  # hasta 3 columnas por fila
columnas = 3

# Crear figura y ejes
fig, axes = plt.subplots(filas, columnas, figsize=(12, filas * 3))
axes = axes.flatten()

# Graficar cada banda
for i, banda in enumerate(bandas):
    ax = axes[i]
    ax.scatter(df[banda], df['ln_turbidez'], alpha=0.6, s=10)
    ax.set_title(f'ln_turbidez vs {banda}', fontsize=10)
    ax.set_xlabel(banda)
    ax.set_ylabel('ln_turbidez')
    ax.grid(True)

# Ocultar los subplots vacíos si hay
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

# Ajustar espacio
plt.tight_layout()
plt.show()

El orden general se mantiene, pero los valores de r bajan un poco al aplicar el logaritmo a la turbidez. Las mejores bandas siguen siendo: B05 > B04 > B06 > B07 > B08

Aunque r haya bajado ligeramente, esas bandas siguen siendo las más prometedoras predictoras lineales.

4.3 Prueba de correlación C3: logaritmo en las bandas

Anteriormente se estaba trabajando con una base de datos de reflectancias procesadas con Sen2Cor. Actualmente trabajamos con una nueva base de datos con las reflectancias corregidas por ACOLITE.

Acolite y Sen2Cor son dos herramientas de software utilizadas para la corrección atmosférica de imágenes satelitales. Acolite se destaca en ambientes acuáticos eutróficos, mientras que Sen2Cor es un procesador más general para la creación de productos Sentinel-2 de nivel 2A.

IMPORTANTE: Este cambio provocó errores en el código a partir de esta sección en adelante, el error se debe a que hay valores de reflectancia que son cero y negativos en la banda 12.

En las pruebas de correlación C3 y C4 aplicamos logaritmo a las bandas. Por defincion de logaritmo no ponemos tomar valores cero ni negativos. Para salvar este error, se procedió a filtrar filas con esos valores.

Creaamos un nuevo DataFrame para aplicarle el logaritmo a las bandas.

Código 
import pandas as pd

# Leer los datos
Datos = pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')

# Guardar cantidad original de filas
filas_originales = Datos.shape[0]

# Filtrar filas: nos quedamos solo con aquellas donde **todos los valores sean mayores a 0**
Datos_filtrado = Datos[(Datos > 0).all(axis=1)]

# Calcular cuántas filas se eliminaron
filas_finales = Datos_filtrado.shape[0]
filas_eliminadas = filas_originales - filas_finales

# Mostrar resultado
print(f"Se eliminaron {filas_eliminadas} filas que contenían ceros o valores negativos.")

# Si querés seguir trabajando con el DataFrame filtrado:
DatosC3_C4 = Datos_filtrado
#DatosC3_C4
Se eliminaron 30 filas que contenían ceros o valores negativos.

Aplicamos logaritmo a las bandas

Código 
import numpy as np

DatosC3 = DatosC3_C4.copy()

col = [col for col in DatosC3.columns if col.startswith('B')]

DatosC3[col] = np.log(DatosC3[col])

#Cambiamos en nombre las columnas, agremamos ln_ a cada columna

DatosC3.columns = ['ln_' + col for col in DatosC3.columns]

DatosC3.rename(columns={'ln_turbidez': 'turbidez'}, inplace=True)

#Creo archivo csv para usarlo para entrenar modelos
DatosC3.to_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C3_turb_vs_ln_banda.csv", index=False)

#DatosC3

Calculamos r entre turb y ln de cada banda, con la función .corr de pandas.

Código 
bandas_ln = [col for col in DatosC3.columns if col.startswith('ln_B')]

correlaciones = {}

for banda in bandas_ln:
    r = DatosC3['turbidez'].corr(DatosC3[banda])
    correlaciones[banda] = r

#Creamos un Data Frame.
df_correlaciones3 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones3 = df_correlaciones3.sort_values(by='r', ascending=False).reset_index(drop=True)

df_correlaciones3['Banda'] = "turb vs " + df_correlaciones3['Banda'].astype(str)
df_correlaciones3.rename(columns={'Banda': 'Combinación 3'}, inplace=True)


df_correlaciones3.to_csv(r'D:\GIT\Turbidez\Datos creados\Correlaciones\C3_turb_vs_ln_banda.csv', index=False)

df_correlaciones3
Combinación 3 r
0 turb vs ln_B05 0.734057
1 turb vs ln_B06 0.717491
2 turb vs ln_B07 0.706154
3 turb vs ln_B04 0.699797
4 turb vs ln_B08 0.698987
5 turb vs ln_B8A 0.641109
6 turb vs ln_B03 0.598708
7 turb vs ln_B02 0.543950
8 turb vs ln_B01 0.533258
9 turb vs ln_B11 0.352344
10 turb vs ln_B12 0.351269

ln_B05 , ln_B06 , ln_B07 , ln_B08

Estas cuatro variables transformadas logarítmicamente se correlacionan de forma fuerte con la turbidez.

Gráfica

Código 
import pandas as pd
import matplotlib.pyplot as plt

# Cargar los datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C3_turb_vs_ln_banda.csv")

# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col.startswith('ln_B')]

# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col.startswith('ln_B')]

# Cantidad de bandas
n = len(bandas)

# Calcular las filas y columnas para el grid de subplots
filas = (n + 2) // 3  # hasta 3 columnas por fila
columnas = 3

# Crear figura y ejes
fig, axes = plt.subplots(filas, columnas, figsize=(12, filas * 3))
axes = axes.flatten()

# Graficar cada banda
for i, banda in enumerate(bandas):
    ax = axes[i]
    ax.scatter(df[banda], df['turbidez'], alpha=0.6, s=10)
    ax.set_title(f'Turbidez vs {banda}', fontsize=10)
    ax.set_xlabel(banda)
    ax.set_ylabel('Turbidez')
    ax.grid(True)

# Ocultar los subplots vacíos si hay
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

# Ajustar espacio
plt.tight_layout()
plt.show()

4.4 Prueba de correlación C4: logaritmo en turbidez y en bandas

Primeramente verificamos que no hayan valores negativos y cero. Ya que por defincion de logaritno no los podremos utilizar.

Código 
import pandas as pd

# Leer los datos
Datos = pd.read_csv(r'D:\GIT\Turbidez\Turbidez_FINAL.csv')

# Guardar cantidad original de filas
filas_originales = Datos.shape[0]

# Filtrar filas: nos quedamos solo con aquellas donde **todos los valores sean mayores a 0**
Datos_filtrado = Datos[(Datos > 0).all(axis=1)]

# Calcular cuántas filas se eliminaron
filas_finales = Datos_filtrado.shape[0]
filas_eliminadas = filas_originales - filas_finales

# Mostrar resultado
print(f"Se eliminaron {filas_eliminadas} filas que contenían ceros o valores negativos.")

# Si querés seguir trabajando con el DataFrame filtrado:
DatosC3_C4 = Datos_filtrado
#DatosC3_C4
Se eliminaron 30 filas que contenían ceros o valores negativos.

Creaamos un nuevo DataFrame para aplicarle el logaritmo a todas las columnas

Importamos nunpy para operar con funciones matemáticas

Código 
import numpy as np

DatosC4 = np.log(DatosC3_C4)
#Cambiamos en nombre las columnas, agremamos ln_ a cada columna
DatosC4.columns = ['ln_' + col for col in DatosC4.columns]

#Creo archivo csv para usarlo para entrenar modelos
DatosC4.to_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv", index=False)

#DatosC4.head()

Calculamos r entre el ln_turb y ln de cada banda, con la función .corr de pandas.

Código 
bandas_ln = [col for col in DatosC4.columns if col.startswith('ln_B')]

correlaciones = {}

for banda in bandas_ln:
    r = DatosC4['ln_turbidez'].corr(DatosC4[banda])
    correlaciones[banda] = r

#Creamos un Data Frame.
df_correlaciones4 = pd.DataFrame(list(correlaciones.items()), columns=['Banda', 'r'])
#Ordenamos por mayor a menos valor de r.
df_correlaciones4 = df_correlaciones4.sort_values(by='r', ascending=False).reset_index(drop=True)

df_correlaciones4['Banda'] = "ln_turb vs " + df_correlaciones4['Banda'].astype(str)
df_correlaciones4.rename(columns={'Banda': 'Combinación 4'}, inplace=True)

df_correlaciones4.to_csv(r'D:\GIT\Turbidez\Datos creados\Correlaciones\C4_ln_turb_vs_ln_banda.csv', index=False)

df_correlaciones4
Combinación 4 r
0 ln_turb vs ln_B05 0.844506
1 ln_turb vs ln_B04 0.787566
2 ln_turb vs ln_B06 0.770889
3 ln_turb vs ln_B07 0.754235
4 ln_turb vs ln_B08 0.740944
5 ln_turb vs ln_B8A 0.661797
6 ln_turb vs ln_B03 0.533858
7 ln_turb vs ln_B02 0.447893
8 ln_turb vs ln_B01 0.397983
9 ln_turb vs ln_B11 0.250618
10 ln_turb vs ln_B12 0.223253

Gráfica

Código 
import pandas as pd
import matplotlib.pyplot as plt

# Cargar los datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")

# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'ln_turbidez']

# Crear una figura con subplots
# Seleccionar las columnas de bandas (excepto turbidez)
bandas = [col for col in df.columns if col != 'ln_turbidez']

# Cantidad de bandas
n = len(bandas)

# Calcular las filas y columnas para el grid de subplots
filas = (n + 2) // 3  # hasta 3 columnas por fila
columnas = 3

# Crear figura y ejes
fig, axes = plt.subplots(filas, columnas, figsize=(12, filas * 3))
axes = axes.flatten()

# Graficar cada banda
for i, banda in enumerate(bandas):
    ax = axes[i]
    ax.scatter(df[banda], df['ln_turbidez'], alpha=0.6, s=10)
    ax.set_title(f'ln_turbidez vs {banda}', fontsize=10)
    ax.set_xlabel(banda)
    ax.set_ylabel('ln_turbidez')
    ax.grid(True)

# Ocultar los subplots vacíos si hay
for j in range(i + 1, len(axes)):
    fig.delaxes(axes[j])

# Ajustar espacio
plt.tight_layout()
plt.show()

B05, B04, B06, B07, B08 son los mejores candidatos como predictores en un modelo de regresión. Tienen una correlación fuerte con la turbidez, por encima de 0.78.

Comparación con combinaciones anteriores:

Combinación Variables transformadas Mejor r
1 turbidez vs Banda B05 (0.88)
2 ln(turb) vs Banda B05 (0.80)
3 turbidez vs ln(Banda) B05 (0.72)
4 ln(turb) vs ln(Banda) B05 (0.87)

5 Seleción de variables

Se resume en una tabla las pruebas de correlación realizadas.

Código 
import pandas as pd 

C1 = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C1_turb_vs_banda.csv")
C2= pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C2_ln_turb_vs_banda.csv")
C3 = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C3_turb_vs_ln_banda.csv")
C4 = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Correlaciones\C4_ln_turb_vs_ln_banda.csv")

C_comb = pd.concat([C1, C2 ,C3, C4], axis=1)

C_comb
Combinación 1 r Combinación 2 r Combinación 3 r Combinación 4 r
0 turb vs B05 0.887299 ln_turb vs B05 0.813367 turb vs ln_B05 0.734057 ln_turb vs ln_B05 0.844506
1 turb vs B06 0.876827 ln_turb vs B04 0.752422 turb vs ln_B06 0.717491 ln_turb vs ln_B04 0.787566
2 turb vs B08 0.869149 ln_turb vs B06 0.714446 turb vs ln_B07 0.706154 ln_turb vs ln_B06 0.770889
3 turb vs B07 0.867741 ln_turb vs B07 0.707941 turb vs ln_B04 0.699797 ln_turb vs ln_B07 0.754235
4 turb vs B04 0.827019 ln_turb vs B08 0.703117 turb vs ln_B08 0.698987 ln_turb vs ln_B08 0.740944
5 turb vs B8A 0.788088 ln_turb vs B8A 0.603662 turb vs ln_B8A 0.641109 ln_turb vs ln_B8A 0.661797
6 turb vs B03 0.673358 ln_turb vs B03 0.514314 turb vs ln_B03 0.598708 ln_turb vs ln_B03 0.533858
7 turb vs B02 0.614315 ln_turb vs B02 0.429268 turb vs ln_B02 0.543950 ln_turb vs ln_B02 0.447893
8 turb vs B01 0.611305 ln_turb vs B01 0.391705 turb vs ln_B01 0.533258 ln_turb vs ln_B01 0.397983
9 turb vs B12 0.429200 ln_turb vs B12 0.242050 turb vs ln_B11 0.352344 ln_turb vs ln_B11 0.250618
10 turb vs B11 0.409447 ln_turb vs B11 0.222510 turb vs ln_B12 0.351269 ln_turb vs ln_B12 0.223253

A partir del análisis con las correlaciones entre turbidez (y su logaritmo) y distintas bandas (y sus logaritmos). Se puede observar que en las combinaciones 1 y 4 las correlaciones son mas altas respecto a C2 y C3. Por que que se seguirá el análisis con estas combinaciones:

• Combinación 1: turb vs bandas

• Combinación 4: ln_turb vs ln_bandas

Se entrenarán modelos con ambas combinaciones por separado y se evalurá su desempeño.

Todos los modelos serán sometidos a bootstrapping con enfoque out-of-bag (OOB):

→ Donde se tomará el 75% de los datos con reemplazo para entrenamiento; → Luego se usa el 25% restante (no muestreados) como conjunto de prueba (test).

Proponemos avanzar hacia un modelo multivariable, utilizando el AIC [11] como criterio de seleccion de variables, este criterio es muy conveniente para encontrar el modelo de mejor balance entre ajuste y complejidad.

¿Qué hace el AIC?

→ Penaliza modelos con demasiadas variables para evitar sobreajuste.

→ Se busca minimizar el AIC: cuanto más bajo, mejor el modelo.

→ No requiere que las variables tengan alta correlación individual, pero sí que mejoren el modelo conjunto.

AIC=n⋅ln(RMSE²)+2k

5.1 Combinación 1: turb vs bandas

Código 
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from utils import run_bootstrap_linear_regression_analysis  # Asegurate de que esta función devuelve lo necesario

# Leer los datos
Datos = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")
y = Datos['turbidez']

# Variables
variables_fijas = ['B05']
variables_a_agregar = ['B06', 'B08', 'B07', 'B04', 'B8A', 'B03', 'B02', 'B01', 'B12', 'B11']

# Resultados
resultados = []

# Configuración de bootstrapping
n_iteraciones_bootstrap = 1000

# Entrenamiento incremental
for i in range(len(variables_a_agregar) + 1):
    variables_usadas = variables_fijas + variables_a_agregar[:i]
    X = Datos[variables_usadas].values

    # División entrenamiento/test
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.25, shuffle=True, random_state=42
    )

    # Bootstrapping: obtener coeficientes promedio y métricas de entrenamiento
    coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = run_bootstrap_linear_regression_analysis(
    X_train, y_train.to_numpy(), n_iteraciones=n_iteraciones_bootstrap)


    # Modelo final con coeficientes promedio
    modelo_final = LinearRegression()
    modelo_final.coef_ = coef_prom
    modelo_final.intercept_ = intercept_prom

    # Predicción sobre testeo
    y_pred = modelo_final.predict(X_test)
    r2 = r2_score(y_test, y_pred)
    rmse = np.sqrt(mean_squared_error(y_test, y_pred))

    # R² ajustado
    n_obs = len(y_test)
    n_vars = X_test.shape[1]
    r2_ajustado = 1 - (1 - r2) * (n_obs - 1) / (n_obs - n_vars - 1)

    # AIC
    residuals = y_test - y_pred
    rss = np.sum(residuals ** 2)
    k = X_test.shape[1] + 1  # +1 por el intercepto
    aic = n_obs * np.log(rss / n_obs) + 2 * k

    # Guardar resultados
    resultados.append({
        "variables": ", ".join(variables_usadas),
        "num_variables": len(variables_usadas),
        "R²_train (bootstrap)": r2_train_boot,
        "RMSE_train (bootstrap)": rmse_train_boot,
        "R²_test": r2,
        "R²-ajustado": r2_ajustado,
        "RMSE_test": rmse,
        "AIC": aic
    })

# Convertir a DataFrame
df_resultados = pd.DataFrame(resultados)
df_resultados
variables num_variables R²_train (bootstrap) RMSE_train (bootstrap) R²_test R²-ajustado RMSE_test AIC
0 B05 1 0.817414 35.380587 0.301995 0.258369 33.294120 130.193709
1 B05, B06 2 0.818655 35.260133 0.333426 0.244549 32.535871 131.364356
2 B05, B06, B08 3 0.825575 34.580912 0.278433 0.123812 33.851391 134.791284
3 B05, B06, B08, B07 4 0.842854 32.823338 0.410136 0.228639 30.606533 133.163686
4 B05, B06, B08, B07, B04 5 0.928892 22.079519 0.502045 0.294564 28.121155 132.114797
5 B05, B06, B08, B07, B04, B8A 6 0.943200 19.733626 0.762204 0.632497 19.433035 120.811081
6 B05, B06, B08, B07, B04, B8A, B03 7 0.947473 18.976730 0.715553 0.516439 21.253927 126.035501
7 B05, B06, B08, B07, B04, B8A, B03, B02 8 0.971103 14.075346 0.860893 0.737242 14.863235 115.160065
8 B05, B06, B08, B07, B04, B8A, B03, B02, B01 9 0.976241 12.762692 0.902274 0.792332 12.457859 110.804661
9 B05, B06, B08, B07, B04, B8A, B03, B02, B01, B12 10 0.979273 11.920504 0.925210 0.818367 10.898344 107.989991
10 B05, B06, B08, B07, B04, B8A, B03, B02, B01, B... 11 0.980290 11.624634 0.917927 0.767460 11.416642 111.662596

Gráfico de AIC

Código 
aic_scores = [r["AIC"] for r in resultados]
n_vars = [r["num_variables"] for r in resultados]

plt.plot(n_vars, aic_scores, marker='o', color='purple')
plt.title('AIC vs Número de Variables')
plt.xlabel('Número de variables')
plt.ylabel('AIC')
plt.grid(True)
plt.show()

¿Por qué el AIC baja al agregar más variables, si se supone que penaliza la complejidad?

AIC=n⋅ln(RMSE²)+2k

donde:

n = número de observaciones k = número de parámetros (incluye el intercepto) ln(RMSE²) representa el ajuste del modelo

AIC no busca modelos simples en sí, sino la mejor compensación entre ajuste y complejidad.

A pesar de agregar más variables, el error (RMSE) está disminuyendo mucho más rápido de lo que penaliza el AIC con +2k

5.2 Combinación 4: ln_turb vs ln_bandas

Código 
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from utils import run_bootstrap_linear_regression_analysis  # Asegurate de que esta función devuelve lo necesario

# Leer los datos
Datos = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")
y = Datos['ln_turbidez']

# Variables
variables_fijas = ['ln_B05']
variables_a_agregar = ['ln_B04', 'ln_B06', 'ln_B07', 'ln_B08', 'ln_B8A', 'ln_B03', 'ln_B02', 'ln_B01','ln_B11' ,'ln_B12']

# Resultados
resultados = []

# Configuración de bootstrapping
n_iteraciones_bootstrap = 1000

# Entrenamiento incremental
for i in range(len(variables_a_agregar) + 1):
    variables_usadas = variables_fijas + variables_a_agregar[:i]
    X = Datos[variables_usadas].values

    # División entrenamiento/test
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.25, shuffle=True, random_state=42
    )

    # Bootstrapping: obtener coeficientes promedio y métricas de entrenamiento
    coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = run_bootstrap_linear_regression_analysis(
    X_train, y_train.to_numpy(), n_iteraciones=n_iteraciones_bootstrap)


    # Modelo final con coeficientes promedio
    modelo_final = LinearRegression()
    modelo_final.coef_ = coef_prom
    modelo_final.intercept_ = intercept_prom

    # Predicción sobre testeo
    y_pred = modelo_final.predict(X_test)
    r2 = r2_score(y_test, y_pred)
    rmse = np.sqrt(mean_squared_error(y_test, y_pred))

    # Cálculo R² ajustado
    n_obs = len(y_test)
    n_vars = X_test.shape[1]

    if n_obs - n_vars - 1 != 0:
     r2_ajustado = 1 - (1 - r2) * (n_obs - 1) / (n_obs - n_vars - 1)
    else:
     r2_ajustado = np.nan  # o podrías usar 0.0 si preferís



    # AIC
    residuals = y_test - y_pred
    rss = np.sum(residuals ** 2)
    k = X_test.shape[1] + 1  # +1 por el intercepto
    aic = n_obs * np.log(rss / n_obs) + 2 * k

    # Guardar resultados
    resultados.append({
        "variables": ", ".join(variables_usadas),
        "num_variables": len(variables_usadas),
        "R²_train (bootstrap)": r2_train_boot,
        "RMSE_train (bootstrap)": rmse_train_boot,
        "R²_test": r2,
        "R²-ajustado": r2_ajustado,
        "RMSE_test": rmse,
        "AIC": aic
    })

# Convertir a DataFrame
df_resultados = pd.DataFrame(resultados)
df_resultados
variables num_variables R²_train (bootstrap) RMSE_train (bootstrap) R²_test R²-ajustado RMSE_test AIC
0 ln_B05 1 0.683149 0.731301 0.772278 0.746976 0.600810 -7.208492
1 ln_B05, ln_B04 2 0.814223 0.559970 0.764594 0.705742 0.610863 -4.843421
2 ln_B05, ln_B04, ln_B06 3 0.911953 0.385503 0.756236 0.651765 0.621613 -2.459633
3 ln_B05, ln_B04, ln_B06, ln_B07 4 0.909160 0.391569 0.732322 0.553871 0.651390 0.569765
4 ln_B05, ln_B04, ln_B06, ln_B07, ln_B08 5 0.908863 0.392208 0.731631 0.463261 0.652231 2.598159
5 ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A 6 0.944083 0.307215 0.559205 -0.101988 0.835898 10.056536
6 ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... 7 0.944316 0.306573 0.579665 -0.401118 0.816268 11.533732
7 ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... 8 0.946340 0.300950 0.639149 -0.804256 0.756309 11.855272
8 ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... 9 0.946290 0.301090 0.727218 -1.727822 0.657571 10.777558
9 ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... 10 0.947751 0.296966 0.710570 NaN 0.677340 13.429195
10 ln_B05, ln_B04, ln_B06, ln_B07, ln_B08, ln_B8A... 11 0.951308 0.286679 0.696386 4.036139 0.693738 15.955473

Gráfico de AIC

Código 
aic_scores = [r["AIC"] for r in resultados]
n_vars = [r["num_variables"] for r in resultados]

plt.plot(n_vars, aic_scores, marker='o', color='purple')
plt.title('AIC vs Número de Variables')
plt.xlabel('Número de variables')
plt.ylabel('AIC')
plt.grid(True)
plt.show()

6 Entrenamiento de Modelos

6.1 Modelos sin logaritmo (turb vs bandas)

6.1.1 Modelo con B05

Código 
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from utils import run_bootstrap_linear_regression_analysis

# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")

# Variables predictoras
variables = ['B05']
X_completo = df[variables].values
y_completo = df['turbidez'].values

# Separar entrenamiento y testeo
X_train, X_test, y_train, y_test = train_test_split(
    X_completo, y_completo, test_size=0.25, random_state=42
)

# Ejecutar bootstrapping
n_iteraciones_config = 1000
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = \
    run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones=n_iteraciones_config)

# ------------------ Guardar coeficientes del modelo B05 ------------------ #
coef_B05 = coef_prom[0]         # coeficiente de la banda B05
intercept_B05 = intercept_prom  # intercepto

# Ecuación en LaTeX (solo para mostrar)
ecuacion_latexB05 = f"turbidez = {coef_B05:.4f} \\cdot B05 + {intercept_B05:.4f}"

from IPython.display import display, Math
display(Math(ecuacion_latexB05))


# Modelo final con coeficientes promedio
modelo_final_promedio = LinearRegression()
modelo_final_promedio.coef_ = coef_prom
modelo_final_promedio.intercept_ = intercept_prom

# Predicción sobre entrenamiento y test
y_train_pred = modelo_final_promedio.predict(X_train)
y_test_pred = modelo_final_promedio.predict(X_test)

# Métricas en testeo
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))

# ------------------ Gráfico con estilo personalizado ------------------ #
plt.figure(figsize=(4.5, 4.5))

# Entrenamiento
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')

# Test
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7, label="Datos de testeo", marker='^')

# Línea ideal
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")

# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")

plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
    f"Regresión lineal con Bootstrapping ({n_iteraciones_config} modelos)\n"
    f"R² entrenamiento: {r2_train_boot:.4f}, RMSE: {rmse_train_boot:.4f} | "
    f"R² testeo: {r2_test:.4f}, RMSE: {rmse_test:.4f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

# ------------------ Coeficientes ------------------ #

# Diccionario para renombrar las variables en notación matemática
nombre_vars_latex = {
    'B05': 'B05'
}

# Construir términos con nombres renombrados
terminos_latex = [
    f"{coef:.4f} \\cdot {nombre_vars_latex[var]}"
    for var, coef in zip(variables, coef_prom)
]

# Construir la ecuación completa en LaTeX
ecuacion_latexB05 = " + ".join(terminos_latex)
ecuacion_latexB05 = f"turbidez = {ecuacion_latexB05} + {intercept_prom:.4f}"

display(Math(ecuacion_latexB05))

\(\displaystyle turbidez = 1286.4244 \cdot B05 + -39.8511\)

\(\displaystyle turbidez = 1286.4244 \cdot B05 + -39.8511\)

6.1.2 Modelo con 6 bandas

Código 
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from utils import run_bootstrap_linear_regression_analysis

# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")

# Variables predictoras
variables = ['B05', 'B06', 'B08', 'B07', 'B04', 'B8A']
X_completo = df[variables].values
y_completo = df['turbidez'].values

# Separar entrenamiento y testeo
X_train, X_test, y_train, y_test = train_test_split(
    X_completo, y_completo, test_size=0.25, random_state=42
)

# Ejecutar bootstrapping
n_iteraciones_config = 1000
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = \
    run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones=n_iteraciones_config)

# Modelo final con coeficientes promedio
modelo_final_promedio = LinearRegression()
modelo_final_promedio.coef_ = coef_prom
modelo_final_promedio.intercept_ = intercept_prom

# Predicción sobre entrenamiento y test
y_train_pred = modelo_final_promedio.predict(X_train)
y_test_pred = modelo_final_promedio.predict(X_test)

# Métricas en testeo
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))

# ------------------ Gráfico con estilo personalizado ------------------ #
plt.figure(figsize=(4.5, 4.5))

# Entrenamiento
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')

# Test
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7, label="Datos de testeo", marker='^')

# Línea ideal
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")

# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")

plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
    f"Regresión lineal con Bootstrapping ({n_iteraciones_config} modelos)\n"
    f"R² entrenamiento: {r2_train_boot:.4f}, RMSE: {rmse_train_boot:.4f} | "
    f"R² testeo: {r2_test:.4f}, RMSE: {rmse_test:.4f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

# ------------------ Coeficientes ------------------ #

# Diccionario para renombrar las variables en notación matemática
nombre_vars_latex = {
    'B05': 'B05',
    'B04': 'B04',
    'B06': 'B06',
    'B07': 'B07',
    'B8A': 'B8A',
    'B08': 'B08'
}

# Construir términos con nombres renombrados
terminos_latex = [
    f"{coef:.4f} \\cdot {nombre_vars_latex[var]}"
    for var, coef in zip(variables, coef_prom)
]

# Construir la ecuación completa en LaTeX
ecuacion_latex = " + ".join(terminos_latex)
ecuacion_latex = f"turbidez = {ecuacion_latex} + {intercept_prom:.4f}"

from IPython.display import display, Math
display(Math(ecuacion_latex))

\(\displaystyle turbidez = 400.9979 \cdot B05 + 5763.2933 \cdot B06 + 1707.9858 \cdot B08 + -1647.7470 \cdot B07 + -1141.7309 \cdot B04 + -3926.9242 \cdot B8A + 2.4209\)

6.2 Modelos aplicando logaritmo en turbidez y en bandas

6.2.1 Modelo con 1 variable (ln_B05)

Escala logarítmica

Código 
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from utils import run_bootstrap_linear_regression_analysis

# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")

# Variables predictoras
variables = ['ln_B05']
X_completo = df[variables].values
y_completo = df['ln_turbidez'].values

# Separar entrenamiento y testeo
X_train, X_test, y_train, y_test = train_test_split(
    X_completo, y_completo, test_size=0.25, random_state=42
)

# Ejecutar bootstrapping
n_iteraciones_config = 1000
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = \
    run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones=n_iteraciones_config)

# Modelo final con coeficientes promedio
modelo_final_promedio = LinearRegression()
modelo_final_promedio.coef_ = coef_prom
modelo_final_promedio.intercept_ = intercept_prom

# Predicción sobre entrenamiento y test
y_train_pred = modelo_final_promedio.predict(X_train)
y_test_pred = modelo_final_promedio.predict(X_test)

# Métricas en testeo
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))

# ------------------ Gráfico con estilo personalizado ------------------ #
plt.figure(figsize=(4.5, 4.5))

# Entrenamiento
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')

# Test
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7, label="Datos de testeo", marker='^')

# Línea ideal
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")

# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")

plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
    f"Regresión lineal con Bootstrapping ({n_iteraciones_config} modelos)\n"
    f"R² entrenamiento: {r2_train_boot:.4f}, RMSE: {rmse_train_boot:.4f} | "
    f"R² testeo: {r2_test:.4f}, RMSE: {rmse_test:.4f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

# ------------------ Coeficientes ------------------ #

# Diccionario para renombrar las variables en notación matemática
nombre_vars_latex = {
    'ln_B05': 'ln(B05)',
    'ln_B04': 'ln(B04)',
    'ln_B06': 'ln(B06)'
}

# Construir términos con nombres renombrados
terminos_latex = [
    f"{coef:.4f} \\cdot {nombre_vars_latex[var]}"
    for var, coef in zip(variables, coef_prom)
]

# Construir la ecuación completa en LaTeX
ecuacion_latex = " + ".join(terminos_latex)
ecuacion_latex = f"\\ln(\\text{{turbidez}}) = {ecuacion_latex} + {intercept_prom:.4f}"

from IPython.display import display, Math
display(Math(ecuacion_latex))

\(\displaystyle \ln(\text{turbidez}) = 1.5714 \cdot ln(B05) + 7.6043\)

Escala real

Código 
# Transformo a escala real

y_train_pred_real = np.exp(y_train_pred)
y_test_pred_real = np.exp(y_test_pred)

y_train_real = np.exp(y_train)
y_test_real = np.exp(y_test)

rmse_test_real = np.sqrt(mean_squared_error(y_test_real, y_test_pred_real))
r2_test_real = r2_score(y_test_real, y_test_pred_real)

#Grafico

plt.figure(figsize=(4.5, 4.5))

# Entrenamiento (en escala real)
plt.scatter(y_train_real, y_train_pred_real, color="#9D50A6", alpha=0.5, label="Entrenamiento")

# Test (en escala real)
plt.scatter(y_test_real, y_test_pred_real, color="red", alpha=0.7, label="Test")

# Línea ideal
min_val = min(y_train_real.min(), y_test_real.min())
max_val = max(y_train_real.max(), y_test_real.max())
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")

# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")

plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
    f"Regresión lineal en escala real\n"
    f"R² testeo: {r2_test_real:.4f}, RMSE: {rmse_test_real:.4f} NTU"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

6.2.2 Modelo con 3 variables

Escala Logarítmica

Código 
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from utils import run_bootstrap_linear_regression_analysis

# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C4_ln_turb_ln_banda.csv")

# Variables predictoras
variables = ['ln_B05','ln_B04','ln_B06']
X_completo = df[variables].values
y_completo = df['ln_turbidez'].values

# Separar entrenamiento y testeo
X_train, X_test, y_train, y_test = train_test_split(
    X_completo, y_completo, test_size=0.25, random_state=42
)

# Ejecutar bootstrapping
n_iteraciones_config = 1000
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = \
    run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones=n_iteraciones_config)

# Modelo final con coeficientes promedio
modelo_final_promedio = LinearRegression()
modelo_final_promedio.coef_ = coef_prom
modelo_final_promedio.intercept_ = intercept_prom

# Predicción sobre entrenamiento y test
y_train_pred = modelo_final_promedio.predict(X_train)
y_test_pred = modelo_final_promedio.predict(X_test)

# Métricas en testeo
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))

# ------------------ Gráfico con estilo personalizado ------------------ #
plt.figure(figsize=(4.5, 4.5))

# Entrenamiento
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')

# Test
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7, label="Datos de testeo", marker='^')

# Línea ideal
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")

# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")

plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
    f"Regresión lineal con Bootstrapping ({n_iteraciones_config} modelos)\n"
    f"R² entrenamiento: {r2_train_boot:.4f}, RMSE: {rmse_train_boot:.4f} | "
    f"R² testeo: {r2_test:.4f}, RMSE: {rmse_test:.4f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

# ------------------ Coeficientes ------------------ #

# Diccionario para renombrar las variables en notación matemática
nombre_vars_latex = {
    'ln_B05': 'ln(B05)',
    'ln_B04': 'ln(B04)',
    'ln_B06': 'ln(B06)'
}

# Construir términos con nombres renombrados
terminos_latex = [
    f"{coef:.4f} \\cdot {nombre_vars_latex[var]}"
    for var, coef in zip(variables, coef_prom)
]

# Construir la ecuación completa en LaTeX
ecuacion_latex = " + ".join(terminos_latex)
ecuacion_latex = f"\\ln(\\text{{turbidez}}) = {ecuacion_latex} + {intercept_prom:.4f}"

from IPython.display import display, Math
display(Math(ecuacion_latex))

# Ecuacion real guardada como string
# Exponenciar el intercepto
intercepto_escala_real = np.exp(intercept_prom)

# Armar términos con potencias
terminos_real = [
    f"{var[3:]}^{{{coef:.4f}}}"  # Elimina "ln_" del nombre de la banda
    for var, coef in zip(variables, coef_prom)
]

# Construir la ecuación en LaTeX como string reutilizable
ecuacion_real_latex = f"\\text{{turbidez}} = {intercepto_escala_real:.4f} \\cdot " + " \\cdot ".join(terminos_real)

# (Opcional) Mostrar
from IPython.display import Math, display
display(Math(ecuacion_real_latex))

# Guardar en otra variable si querés reutilizar luego como string sin formato LaTeX
ecuacion_turb = f"turbidez = {intercepto_escala_real:.4f} * " + " * ".join([
    f"{var[3:]}^{coef:.4f}" for var, coef in zip(variables, coef_prom)
])

\(\displaystyle \ln(\text{turbidez}) = 8.4385 \cdot ln(B05) + -5.4953 \cdot ln(B04) + -1.8635 \cdot ln(B06) + 6.1811\)

\(\displaystyle \text{turbidez} = 483.5190 \cdot B05^{8.4385} \cdot B04^{-5.4953} \cdot B06^{-1.8635}\)

Escala Real

Código 
# Transformo a escala real

y_train_pred_real = np.exp(y_train_pred)
y_test_pred_real = np.exp(y_test_pred)

y_train_real = np.exp(y_train)
y_test_real = np.exp(y_test)

rmse_test_real = np.sqrt(mean_squared_error(y_test_real, y_test_pred_real))
r2_test_real = r2_score(y_test_real, y_test_pred_real)

#Grafico

plt.figure(figsize=(4.5, 4.5))

# Entrenamiento (en escala real)
plt.scatter(y_train_real, y_train_pred_real, color="#9D50A6", alpha=0.5, label="Entrenamiento")

# Test (en escala real)
plt.scatter(y_test_real, y_test_pred_real, color="red", alpha=0.7, label="Test")

# Línea ideal
min_val = min(y_train_real.min(), y_test_real.min())
max_val = max(y_train_real.max(), y_test_real.max())
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")

# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")

plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
    f"Regresión lineal en escala real\n"
    f"R² testeo: {r2_test_real:.4f}, RMSE: {rmse_test_real:.4f} NTU"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

7 Procesamiento de productos satelitales

Para el procesamiento se utilizará la librería rasterio en python.

Tabla de referencia para seleccionar bandas en rasterio

Nombre de banda en archivo .tif Banda de Sentinel-2 Número para la lectura en rasterio
Banda 1 B01 1
Banda 2 B02 2
Banda 3 B03 3
Banda 4 B04 4
Banda 5 B05 5
Banda 6 B06 6
Banda 7 B07 7
Banda 8 B08 8
Banda 9 B8A 9
Banda 10 B11 10
Banda 11 B12 11

Lectura de imágen raster y visualización del área de estudio, con el objetivo de verificar que el código está funcionando y de que se leen correctamente los archivos.

Código 
import rasterio as rs
import matplotlib.pyplot as plt
import numpy as np

ruta = r"D:\GIT\Turbidez\recorte_acolite\2025-06-09.tif"

# Abrir el archivo
with rs.open(ruta) as raster:

    #Definimos las bandas 
    #Bandas RGB para ver área de estudio
    B02 = raster.read(2) #Azul
    B03 = raster.read(3) #Verde
    B04 = raster.read(4) #Rojo

    #Para calcular NDWI con B03 y B11
    B11 = raster.read(10)

    #para el modelo
    B05 = raster.read(5)
    
# Composión de imagen 
rgb = np.stack([B04, B03, B02], axis=-1)

# Mejorar constraste
def stretch(banda):
    p2, p98 = np.percentile(banda, (2, 98))
    return np.clip((banda - p2) / (p98 - p2 + 1e-6), 0, 1)

rgb_stretched = np.stack([stretch(B04), stretch(B03), stretch(B02)], axis=-1)

# Gráfico
plt.figure(figsize=(4.5, 4.5))
plt.imshow(rgb_stretched)
plt.title("Río Paraná")
plt.axis("off")
plt.tight_layout()
plt.show()

Se utilizrá el índice NDWI para recortar la imágen satelital, y poder trabajar únicamente con los píxeles de agua.

\(NDWI=(B03-B11)/(B03+B11)\)

Código 
import numpy as np
import matplotlib.pyplot as plt
from skimage.filters import threshold_otsu

# NDWI
ndwi = (B03 - B11) / (B03 + B11 + 1e-6)

# Calcular umbral automático con Otsu
umbral = threshold_otsu(ndwi[np.isfinite(ndwi)])  # descartamos NaN e inf

# Binarización: agua = 1, no agua = 0
ndwi_bin = np.where(ndwi >= umbral, 1, 0)

# Gráfico máscara de agua
fig, ax = plt.subplots(figsize=(4.5, 4.5))
img = ax.imshow(ndwi_bin, cmap='Greys')
ax.set_title(f"NDWI - Otsu (umbral={umbral:.3f})")
cbar = fig.colorbar(img, ax=ax, fraction=0.046, pad=0.04)
cbar.set_label("Agua (1) / No Agua (0)")
ax.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False)
ax.axis("on")
plt.tight_layout()
plt.show()

Luego del anális de las distantas pruebas de correlación con diferenctes bandas y combinaciones de ellas, optamos por tomar el modelo que utiliza la banda B05, para estimar la turbidez.

Aplicamos la ecuación hallada para ver la distribución de la turbidez en el río.

Código 
import numpy as np
import matplotlib.pyplot as plt
import re

# --- Calcular turbidez directamente con variables guardadas ---
turbidez = coef_B05 * B05 + intercept_B05

# Forzar a que no existan valores negativos (los pone en 0)
turbidez = np.clip(turbidez, 0, None)

# --- Aplicar máscara NDWI (solo donde hay agua) ---
turbidez_agua = np.where(ndwi_bin == 1, turbidez, np.nan)

# --- Calcular percentiles para ajustar visualización ---
p5, p95 = np.nanpercentile(turbidez_agua, [5, 95])

Mapa de turbidez

Código 
import matplotlib.pyplot as plt

# Crear figura
fig, ax = plt.subplots(figsize=(4.5, 4.5))

# Mostrar imagen con coordenadas reales
img = ax.imshow(turbidez_agua, cmap="rainbow", vmin=p5, vmax=p95)

# Título y colorbar
ax.set_title("Turbidez 2025-06-09")
cbar = fig.colorbar(img, ax=ax, fraction=0.046, pad=0.04)
cbar.set_label("Turbidez (NTU)")
ax.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False)
ax.axis("on")
plt.tight_layout()
plt.show()

Para exportar tif

Código 
import rasterio
import numpy as np
import os
import re

# Ruta de una banda original para copiarle la georreferenciación
ruta_base = r"D:\GIT\Turbidez\recorte_acolite\2025-06-09.tif"

# Carpeta donde guardar el resultado
carpeta_salida = r"D:\GIT\Turbidez\tif salidas"

# --- Extraer la fecha del nombre del archivo base ---
nombre_base = os.path.basename(ruta_base)  # "2025-06-09.tif"
fecha = re.search(r"\d{4}-\d{2}-\d{2}", nombre_base).group()  # "2025-06-09"

# Nombre final con fecha incluida
nombre_salida = f"turbidez_{fecha}.tif"
ruta_salida = os.path.join(carpeta_salida, nombre_salida)

# Abrimos la banda base para obtener el perfil (tamaño, proyección, etc.)
with rasterio.open(ruta_base) as src:
    perfil = src.profile

# Ajustar perfil para nuestro resultado
perfil.update(
    dtype=rasterio.float32,
    count=1,
    compress="lzw",
    nodata=np.nan
)

# Guardar turbidez_agua en un GeoTIFF
with rasterio.open(ruta_salida, "w", **perfil) as dst:
    dst.write(turbidez_agua.astype(rasterio.float32), 1)

8 Lectura automática de tif y procesamiento para mapas de turbidez, cálculo de umbral de mácara de agua con método Otsu

Código 
import rasterio as rs
import numpy as np
import matplotlib.pyplot as plt
import glob
import os
from skimage.filters import threshold_otsu


# Carpeta con archivos .tif
carpeta = r"D:\GIT\Turbidez\recorte_acolite"

# Buscar todos los .tif
tifs = sorted(glob.glob(os.path.join(carpeta, "*.tif")))
print(f"Se encontraron {len(tifs)} archivos .tif")

for ruta in tifs:
    nombre = os.path.basename(ruta)
    fecha = os.path.splitext(nombre)[0]   # nombre sin extensión
    print(f"\nProcesando: {nombre}")

    with rs.open(ruta) as raster:
        # Bandas necesarias
        B03 = raster.read(3).astype(np.float32)  # Verde
        B05 = raster.read(5).astype(np.float32)  # RedEdge (para ecuación turbidez)
        B11 = raster.read(10).astype(np.float32) # SWIR

    # --- NDWI ---
    ndwi = (B03 - B11) / (B03 + B11 + 1e-6)

    # --- Umbral Otsu ---
    umbral = threshold_otsu(ndwi[np.isfinite(ndwi)])
    print(f"   Umbral NDWI (Otsu): {umbral:.3f}")

    # --- Binarización: agua = 1, no agua = 0 ---
    ndwi_bin = np.where(ndwi >= umbral, 1, 0)

    # --- Turbidez ---
    turbidez = coef_B05 * B05 + intercept_B05
    turbidez = np.clip(turbidez, 0, None)  # negativos a 0
    turbidez_agua = np.where(ndwi_bin == 1, turbidez, np.nan)

    # --- Calcular percentiles para ajustar visualización ---
    p5, p95 = np.nanpercentile(turbidez_agua, [5, 95])

    # --- Gráfico conjunto ---
    fig, axes = plt.subplots(1, 2, figsize=(8, 4.5))

    # Panel 1: máscara de agua
    im1 = axes[0].imshow(ndwi_bin, cmap="Greys")
    axes[0].set_title(f"Máscara de agua - {fecha}")
    axes[0].axis("off")
    cbar1 = fig.colorbar(im1, ax=axes[0], fraction=0.046, pad=0.04)
    cbar1.set_label("Agua (1) / No Agua (0)")
    

    # Panel 2: turbidez
    im2 = axes[1].imshow(turbidez_agua, cmap="rainbow", vmin=p5, vmax=p95)
    axes[1].set_title(f"Turbidez - {fecha}")
    axes[1].axis("off")
    cbar2 = fig.colorbar(im2, ax=axes[1], fraction=0.046, pad=0.04)
    cbar2.set_label("Turbidez (NTU)")


    plt.tight_layout()
    plt.show()
Se encontraron 7 archivos .tif

Procesando: 2023-05-11.tif
   Umbral NDWI (Otsu): 0.308


Procesando: 2023-12-12.tif
   Umbral NDWI (Otsu): -0.084


Procesando: 2024-05-20.tif
   Umbral NDWI (Otsu): 37.771


Procesando: 2024-05-30.tif
   Umbral NDWI (Otsu): -0.148


Procesando: 2024-12-16.tif
   Umbral NDWI (Otsu): -138.005


Procesando: 2025-06-09.tif
   Umbral NDWI (Otsu): 0.279


Procesando: 2025-09-12.tif
   Umbral NDWI (Otsu): 1.153

Como el método Otsu no puede procesar correctamente algunas imágenes, se decidió utilizar el valor promedio del NDWI para calcular el umbral.

9 Conclusión para modelos entrenados con métodos de regresion lineal clásicos

Se seleccionará el mejor modelo lineal, este es el que utiliza la B05 de Sentinel-2 como variable predictora de la turbidez y será empleado para hacer la comparación con los métodos de aprendizaje automático próximamente. A continuación se resume la ecuación, métricas de desempeño y mapas de turbidez:

Modelo con B05 y sus métricas

Código 
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from utils import run_bootstrap_linear_regression_analysis

# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")

# Variables predictoras
variables = ['B05']
X_completo = df[variables].values
y_completo = df['turbidez'].values

# Separar entrenamiento y testeo
X_train, X_test, y_train, y_test = train_test_split(
    X_completo, y_completo, test_size=0.25, random_state=42
)

# Ejecutar bootstrapping
n_iteraciones_config = 1000
coef_prom, intercept_prom, r2_train_boot, rmse_train_boot = \
    run_bootstrap_linear_regression_analysis(X_train, y_train, n_iteraciones=n_iteraciones_config)

# ------------------ Guardar coeficientes del modelo B05 ------------------ #
coef_B05 = coef_prom[0]         # coeficiente de la banda B05
intercept_B05 = intercept_prom  # intercepto

# Ecuación en LaTeX (solo para mostrar)
ecuacion_latexB05 = f"turbidez = {coef_B05:.4f} \\cdot B05 + {intercept_B05:.4f}"

from IPython.display import display, Math
display(Math(ecuacion_latexB05))


# Modelo final con coeficientes promedio
modelo_final_promedio = LinearRegression()
modelo_final_promedio.coef_ = coef_prom
modelo_final_promedio.intercept_ = intercept_prom

# Predicción sobre entrenamiento y test
y_train_pred = modelo_final_promedio.predict(X_train)
y_test_pred = modelo_final_promedio.predict(X_test)

# Métricas en testeo
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))

# ------------------ Gráfico con estilo personalizado ------------------ #
plt.figure(figsize=(4.5, 4.5))

# Entrenamiento
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5, label="Datos de entrenamiento", marker='o')

# Test
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7, label="Datos de testeo", marker='^')

# Línea ideal
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--', color="#17A77E", lw=2, label="Línea ideal")

# Línea de tendencia sobre test
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # [pendiente, intercepto]
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black', lw=2, label="Línea de tendencia")

plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")
plt.title(
    f"Regresión lineal con Bootstrapping ({n_iteraciones_config} modelos)\n"
    f"R² entrenamiento: {r2_train_boot:.4f}, RMSE: {rmse_train_boot:.4f} | "
    f"R² testeo: {r2_test:.4f}, RMSE: {rmse_test:.4f}"
)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

# ------------------ Coeficientes ------------------ #

# Diccionario para renombrar las variables en notación matemática
nombre_vars_latex = {
    'B05': 'B05'
}

# Construir términos con nombres renombrados
terminos_latex = [
    f"{coef:.4f} \\cdot {nombre_vars_latex[var]}"
    for var, coef in zip(variables, coef_prom)
]

# Construir la ecuación completa en LaTeX
ecuacion_latexB05 = " + ".join(terminos_latex)
ecuacion_latexB05 = f"turbidez = {ecuacion_latexB05} + {intercept_prom:.4f}"

display(Math(ecuacion_latexB05))

\(\displaystyle turbidez = 1294.7081 \cdot B05 + -40.1019\)

\(\displaystyle turbidez = 1294.7081 \cdot B05 + -40.1019\)

Mapas de turbidez

Código 
import rasterio as rs
import numpy as np
import matplotlib.pyplot as plt
import glob
import os
from skimage.filters import threshold_otsu

# Carpeta con archivos .tif
carpeta = r"D:\GIT\Turbidez\recorte_acolite"
salida = os.path.join(carpeta, "turbidez_out")  # subcarpeta para salidas

# Crear carpeta si no existe
os.makedirs(salida, exist_ok=True)

# Buscar todos los .tif
tifs = sorted(glob.glob(os.path.join(carpeta, "*.tif")))
print(f"Se encontraron {len(tifs)} archivos .tif")

for ruta in tifs:
    nombre = os.path.basename(ruta)
    fecha = os.path.splitext(nombre)[0]   # nombre sin extensión
    print(f"\nProcesando: {nombre}")

    with rs.open(ruta) as raster:
        # Bandas necesarias
        B03 = raster.read(3).astype(np.float32)  # Verde
        B05 = raster.read(5).astype(np.float32)  # RedEdge
        B11 = raster.read(10).astype(np.float32) # SWIR

        # Guardar perfil original
        perfil = raster.profile

    # --- NDWI ---
    ndwi = (B03 - B11) / (B03 + B11 + 1e-6)

    # Calcular el promedio (ignorando NaN si hay)
    umbral = np.nanmean(ndwi)

    # Binarizar con la media
    mascara = (ndwi > umbral).astype(np.uint8)

    # --- Turbidez ---
    turbidez = coef_B05 * B05 + intercept_B05
    turbidez = np.clip(turbidez, 0, None)  # negativos a 0
    turbidez_agua = np.where(mascara == 1, turbidez, np.nan)

    # --- Exportar a GeoTIFF ---
    salida_tif = os.path.join(salida, f"turbidez_{fecha}.tif")

    perfil.update(
        dtype=rs.float32,
        count=1,  # solo una banda (turbidez)
        compress='lzw',
        nodata=np.nan
    )

    with rs.open(salida_tif, 'w', **perfil) as dst:
        dst.write(turbidez_agua.astype(np.float32), 1)


    # --- Calcular percentiles para ajustar visualización ---
    p5, p95 = np.nanpercentile(turbidez_agua, [5, 95])

    # --- Gráfico conjunto ---
    fig, axes = plt.subplots(1, 2, figsize=(8, 4.5))

    # Panel 1: máscara de agua
    im1 = axes[0].imshow(mascara, cmap="Greys")
    axes[0].set_title(f"Máscara de agua - {fecha}")
    axes[0].axis("off")
    cbar1 = fig.colorbar(im1, ax=axes[0], fraction=0.046, pad=0.04)
    cbar1.set_label("Agua (1) / No Agua (0)")

    # Panel 2: turbidez
    im2 = axes[1].imshow(turbidez_agua, cmap="rainbow", vmin=p5, vmax=p95)
    axes[1].set_title(f"Turbidez - {fecha}")
    axes[1].axis("off")
    cbar2 = fig.colorbar(im2, ax=axes[1], fraction=0.046, pad=0.04)
    cbar2.set_label("Turbidez (NTU)")

    plt.tight_layout()
    plt.show()
Se encontraron 7 archivos .tif

Procesando: 2023-05-11.tif


Procesando: 2023-12-12.tif


Procesando: 2024-05-20.tif


Procesando: 2024-05-30.tif


Procesando: 2024-12-16.tif


Procesando: 2025-06-09.tif


Procesando: 2025-09-12.tif

10 Método de aprendizaje automático XGBoost

Código 
from xgboost import XGBRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error
import pandas as pd
import numpy as np

# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")

# Definimos variables predictoras (bandas) y objetivo (turbidez medida)
X = df[["B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08", "B11", "B12", "B8A"]]
y = df["turbidez"]

# Dividir en train/test
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.25, random_state=42)

# Modelo XGBoost
model = XGBRegressor(
    n_estimators=500,
    learning_rate=0.05,
    max_depth=6,
    subsample=0.8,
    colsample_bytree=0.8,
    random_state=42)

# Entrenar
model.fit(X_train, y_train)

# --- Predicciones ---
y_train_pred = model.predict(X_train)
y_test_pred = model.predict(X_test)

# --- Métricas entrenamiento ---
r2_train = r2_score(y_train, y_train_pred)
rmse_train = np.sqrt(mean_squared_error(y_train, y_train_pred))

# --- Métricas testeo ---
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))

# --- Print resultados ---
print("=== Resultados Entrenamiento ===")
print(f"R² train: {r2_train:.3f}")
print(f"RMSE train: {rmse_train:.3f}")

print("\n=== Resultados Testeo ===")
print(f"R² test: {r2_test:.3f}")
print(f"RMSE test: {rmse_test:.3f}")

# Gráfico estimado vs real

import matplotlib.pyplot as plt
import numpy as np

plt.figure(figsize=(4.5, 4.5))

# --- Puntos entrenamiento ---
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5,
            label="Datos de entrenamiento", marker='o')

# --- Puntos test ---
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7,
            label="Datos de testeo", marker='^')

# --- Línea ideal 1:1 ---
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--',
         color="#17A77E", lw=2, label="Línea ideal")

# --- Línea de tendencia sobre test ---
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # pendiente, intercepto
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black',
         lw=2, label="Línea de tendencia (test)")

# --- Etiquetas ---
plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")

plt.title(
    f"Modelo XGBoost\n"
    f"R² train: {r2_train:.3f}, RMSE: {rmse_train:.3f} | "
    f"R² test: {r2_test:.3f}, RMSE: {rmse_test:.3f}"
)

plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()
=== Resultados Entrenamiento ===
R² train: 1.000
RMSE train: 0.007

=== Resultados Testeo ===
R² test: 0.906
RMSE test: 12.208

10.1 Mapas con XGBoost

Código 
import rasterio as rs
import numpy as np
import matplotlib.pyplot as plt
import glob
import os
from skimage.filters import threshold_otsu
import pandas as pd

# --- Configuración ---
carpeta = r"D:\GIT\Turbidez\recorte_acolite"
salida = os.path.join(carpeta, "turbidez_out_XBOOST")
os.makedirs(salida, exist_ok=True)

# Definir nombres de features (idénticos a entrenamiento)
feature_names = ["B01", "B02", "B03", "B04", "B05",
                 "B06", "B07", "B08", "B11", "B12", "B8A"]

# Buscar todos los .tif
tifs = sorted(glob.glob(os.path.join(carpeta, "*.tif")))
print(f"Se encontraron {len(tifs)} archivos .tif")

for ruta in tifs:
    nombre = os.path.basename(ruta)
    fecha = os.path.splitext(nombre)[0]
    print(f"\nProcesando: {nombre}")

    with rs.open(ruta) as raster:
        # Leer bandas
        B01 = raster.read(1).astype(np.float32)
        B02 = raster.read(2).astype(np.float32)
        B03 = raster.read(3).astype(np.float32)
        B04 = raster.read(4).astype(np.float32)
        B05 = raster.read(5).astype(np.float32)
        B06 = raster.read(6).astype(np.float32)
        B07 = raster.read(7).astype(np.float32)
        B08 = raster.read(8).astype(np.float32)
        B8A = raster.read(9).astype(np.float32)
        B11 = raster.read(10).astype(np.float32)
        B12 = raster.read(11).astype(np.float32)

        perfil = raster.profile

    # --- NDWI y máscara de agua ---
    ndwi = (B03 - B11) / (B03 + B11 + 1e-6)
    umbral = np.nanmean(ndwi)
    mascara = (ndwi > umbral).astype(np.uint8)

    # --- Crear array de features ---
    filas, cols = B03.shape
    X_img = np.stack([B01, B02, B03, B04, B05,
                      B06, B07, B08, B11, B12, B8A], axis=-1)
    X_img_2d = X_img.reshape(-1, X_img.shape[-1])

    # Convertir a DataFrame con nombres de columnas iguales al entrenamiento
    X_img_2d_df = pd.DataFrame(X_img_2d, columns=feature_names)

    # --- Predecir turbidez ---
    turbidez_pred = model.predict(X_img_2d_df)
    turbidez_pred = turbidez_pred.reshape(filas, cols)

    # Aplicar máscara de agua
    turbidez_agua = np.where(mascara == 1, turbidez_pred, np.nan)

    # --- Exportar a GeoTIFF ---
    salida_tif = os.path.join(salida, f"turbidez_{fecha}.tif")
    perfil.update(dtype=rs.float32, count=1, compress='lzw', nodata=np.nan)
    with rs.open(salida_tif, 'w', **perfil) as dst:
        dst.write(turbidez_agua.astype(np.float32), 1)

    # --- Percentiles para visualización ---
    p5, p95 = np.nanpercentile(turbidez_agua, [5, 95])

    # --- Gráfico ---
    fig, axes = plt.subplots(1, 2, figsize=(8, 4.5))
    im1 = axes[0].imshow(mascara, cmap="Greys")
    axes[0].set_title(f"Máscara de agua - {fecha}")
    axes[0].axis("off")
    cbar1 = fig.colorbar(im1, ax=axes[0], fraction=0.046, pad=0.04)
    cbar1.set_label("Agua (1) / No Agua (0)")

    im2 = axes[1].imshow(turbidez_agua, cmap="rainbow", vmin=p5, vmax=p95)
    axes[1].set_title(f"Turbidez (XGBoost) - {fecha}")
    axes[1].axis("off")
    cbar2 = fig.colorbar(im2, ax=axes[1], fraction=0.046, pad=0.04)
    cbar2.set_label("Turbidez (NTU)")

    plt.tight_layout()
    plt.show()
Se encontraron 7 archivos .tif

Procesando: 2023-05-11.tif


Procesando: 2023-12-12.tif


Procesando: 2024-05-20.tif


Procesando: 2024-05-30.tif


Procesando: 2024-12-16.tif


Procesando: 2025-06-09.tif


Procesando: 2025-09-12.tif

10.2 Hiperparámetros utilizados

• n_estimators=500 → número de arboles que se entrenan secuencialmente

• learning_rate=0.05 → indica que tanta aprende cada arbol

Es el “paso” con que el modelo corrige los errores. -Valores bajos (0.01–0.1) → aprendizaje más lento pero más estable. -Valores altos (>0.3) → puede sobreajustar rápido.

• max_depth=6 → Profundidad máxima de cada árbol.

Controla la complejidad de las interacciones entre las variables (bandas). Más profundo → aprende patrones complejos, pero también ruido.

• subsample=0.8 → Fracción de los datos usados para entrenar cada árbol.

Introduce aleatoriedad. Ayuda a reducir el sobreajuste, 0.8 significa que cada árbol usa el 80% de las filas del dataset.

• colsample_bytree=0.8 → Fracción de las variables (columnas) usadas para cada árbol. Igual que subsample, pero sobre las bandas.

Esto fuerza a los árboles a usar combinaciones distintas de bandas, lo que mejora la generalización.

0.8 = buena práctica para evitar que el modelo dependa demasiado de una sola banda.

Los hiperparámetros controlan el comportamiento del modelo, pero no se aprenden automáticamente del conjunto de datos. Por eso, necesitamos un método que busque las mejores combinaciones posibles para obtener el máximo rendimiento sin sobreajustar.

10.3 Optimización de hiperparámetros

Existen muchos métodos de optimizacion de hiperparámetros

  1. Grid Search (búsqueda en rejilla)

Prueba todas las combinaciones posibles de los valores que uno defina.

Por ejemplo, si se prueba 3 valores de max_depth y 3 de learning_rate, va a entrenar 3 × 3 = 9 modelos.

Usa validación cruzada (cross-validation) para evaluar cada combinación.

Ventajas

• Simple, exhaustivo y fácil de entender.

• Ideal cuando el número de combinaciones es pequeño.

Desventajas

• Muy lento si probás muchos parámetros o valores.

• No aprende de los resultados anteriores (prueba todo a ciegas).

  1. Random Search (búsqueda aleatoria)

En lugar de probar todas las combinaciones, elige al azar un número fijo de combinaciones.

Ejemplo: podrías probar 30 combinaciones elegidas aleatoriamente de un rango mucho más grande.

Ventajas

• Mucho más rápido que Grid Search.

• Puede explorar un espacio grande de parámetros.

• Buen punto medio entre precisión y costo computacional.

Desventajas

• Resultados algo aleatorios (aunque reproducibles con random_state).

• No garantiza probar la “mejor” combinación posible.

  1. Búsqueda bayesiana inteligente (Optuna / Hyperopt / scikit-optimize)

Usa inteligencia adaptativa: aprende de las pruebas anteriores y dirige las siguientes pruebas hacia las regiones más prometedoras.

Es una búsqueda inteligente, no aleatoria ni exhaustiva.

Ventajas

• Mucho más eficiente que Grid o Random Search.

• Encuentra buenos resultados con menos intentos.

• Ideal para modelos complejos como XGBoost.

Desventajas

• Requiere instalar librerías externas (como optuna).

• Más técnica, aunque muy potente.

10.4 Modelo XGBoost optimizado con Grid Search (Búsqueda por cuadrícula)

Código 
from xgboost import XGBRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error
import pandas as pd
import numpy as np

# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")

# Definimos variables predictoras (bandas) y objetivo (turbidez medida)
X = df[["B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08", "B11", "B12", "B8A"]]
y = df["turbidez"]

# Dividir en train/test
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.25, random_state=42)

# --- OPTIMIZACIÓN DE HIPERPARÁMETROS CON GRID SEARCH ---

from sklearn.model_selection import GridSearchCV
from xgboost import XGBRegressor

# Definimos las combinaciones a probar
param_grid = {
    'max_depth': [4, 6, 8],
    'learning_rate': [0.01, 0.05, 0.1],
    'subsample': [0.7, 0.8, 1.0],
    'colsample_bytree': [0.7, 0.8, 1.0],
    'n_estimators': [300, 500, 700]
}

# Modelo base
xgb_model = XGBRegressor(random_state=42)

# Configuración del GridSearch
grid_search = GridSearchCV(
    estimator=xgb_model,
    param_grid=param_grid,
    scoring='r2',
    cv=5,           # validación cruzada con 5 particiones
    verbose=0,
    n_jobs=1       )

# Entrenamos el GridSearch

grid_search.fit(X_train, y_train)

# Mostramos los mejores parámetros y resultados
print("\n=== RESULTADOS GRID SEARCH ===")
print("Mejores parámetros encontrados:")
import pandas as pd
from IPython.display import display

# Convertir los mejores parámetros en DataFrame
best_params_df = pd.DataFrame([grid_search.best_params_])

# Mostrar tabla
display(best_params_df)

# Guardamos el mejor modelo
model = grid_search.best_estimator_

# Entrenamos nuevamente con los mejores parámetros (opcional)
model.fit(X_train, y_train)

# --- Predicciones ---
y_train_pred = model.predict(X_train)
y_test_pred = model.predict(X_test)

# --- Métricas entrenamiento ---
r2_train = r2_score(y_train, y_train_pred)
rmse_train = np.sqrt(mean_squared_error(y_train, y_train_pred))

# --- Métricas testeo ---
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))

# --- Print resultados ---
print("=== Resultados Entrenamiento ===")
print(f"R² train: {r2_train:.3f}")
print(f"RMSE train: {rmse_train:.3f}")

print("\n=== Resultados Testeo ===")
print(f"R² test: {r2_test:.3f}")
print(f"RMSE test: {rmse_test:.3f}")

# Gráfico estimado vs real

import matplotlib.pyplot as plt
import numpy as np

plt.figure(figsize=(4.5, 4.5))

# --- Puntos entrenamiento ---
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5,
            label="Datos de entrenamiento", marker='o')

# --- Puntos test ---
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7,
            label="Datos de testeo", marker='^')

# --- Línea ideal 1:1 ---
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--',
         color="#17A77E", lw=2, label="Línea ideal")

# --- Línea de tendencia sobre test ---
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # pendiente, intercepto
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black',
         lw=2, label="Línea de tendencia (test)")

# --- Etiquetas ---
plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")

plt.title(
    f"Modelo XGBoost (Optimizado con Grid Searh)\n"
    f"R² train: {r2_train:.3f}, RMSE: {rmse_train:.3f} | "
    f"R² test: {r2_test:.3f}, RMSE: {rmse_test:.3f}"
)

plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

=== RESULTADOS GRID SEARCH ===
Mejores parámetros encontrados:
colsample_bytree learning_rate max_depth n_estimators subsample
0 0.7 0.05 4 300 0.7 
=== Resultados Entrenamiento ===
R² train: 1.000
RMSE train: 0.595

=== Resultados Testeo ===
R² test: 0.904
RMSE test: 12.366

10.5 Modelo XGBoost optimizado con Optuna (Búsqueda bayesiana inteligente)

Código 
from xgboost import XGBRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error
import pandas as pd
import numpy as np

# Leer datos
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")

# Definimos variables predictoras (bandas) y objetivo (turbidez medida)
X = df[["B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08", "B11", "B12", "B8A"]]
y = df["turbidez"]

# Dividir en train/test
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.25, random_state=42)

# --- OPTIMIZACIÓN DE HIPERPARÁMETROS CON OPTUNA ---

import optuna
optuna.logging.set_verbosity(optuna.logging.WARNING)
from xgboost import XGBRegressor
from sklearn.model_selection import cross_val_score

def objective(trial):
    # Espacio de búsqueda de hiperparámetros
    params = {
        'max_depth': trial.suggest_int('max_depth', 3, 10),
        'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.3),
        'subsample': trial.suggest_float('subsample', 0.6, 1.0),
        'colsample_bytree': trial.suggest_float('colsample_bytree', 0.6, 1.0),
        'n_estimators': trial.suggest_int('n_estimators', 200, 800),
        'random_state': 42
    }

    # Modelo con los parámetros actuales
    model = XGBRegressor(**params)

    # Validación cruzada de 5 particiones
    score = cross_val_score(model, X_train, y_train, scoring='r2', cv=5).mean()
    return score

# Crear estudio de optimización
study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=50, show_progress_bar=True)

# Mostrar resultados
print("\n=== RESULTADOS OPTUNA ===")
print("Mejores parámetros encontrados:")
import pandas as pd
from IPython.display import display

# Convertir los mejores parámetros a DataFrame
best_params_df = pd.DataFrame([study.best_params])

# Agregar la métrica objetivo (R²)
best_params_df["R2_promedio_CV"] = study.best_value

# Mostrar tabla
display(best_params_df)





# Entrenar el modelo final con los mejores parámetros
model = XGBRegressor(**study.best_params)
model.fit(X_train, y_train)


# --- Predicciones ---
y_train_pred = model.predict(X_train)
y_test_pred = model.predict(X_test)

# --- Métricas entrenamiento ---
r2_train = r2_score(y_train, y_train_pred)
rmse_train = np.sqrt(mean_squared_error(y_train, y_train_pred))

# --- Métricas testeo ---
r2_test = r2_score(y_test, y_test_pred)
rmse_test = np.sqrt(mean_squared_error(y_test, y_test_pred))

# --- Print resultados ---
print("=== Resultados Entrenamiento ===")
print(f"R² train: {r2_train:.3f}")
print(f"RMSE train: {rmse_train:.3f}")

print("\n=== Resultados Testeo ===")
print(f"R² test: {r2_test:.3f}")
print(f"RMSE test: {rmse_test:.3f}")

# Gráfico estimado vs real

import matplotlib.pyplot as plt
import numpy as np

plt.figure(figsize=(4.5, 4.5))

# --- Puntos entrenamiento ---
plt.scatter(y_train, y_train_pred, color="#9D50A6", alpha=0.5,
            label="Datos de entrenamiento", marker='o')

# --- Puntos test ---
plt.scatter(y_test, y_test_pred, color="red", alpha=0.7,
            label="Datos de testeo", marker='^')

# --- Línea ideal 1:1 ---
min_val = min(np.min(y_train), np.min(y_test))
max_val = max(np.max(y_train), np.max(y_test))
plt.plot([min_val, max_val], [min_val, max_val], '--',
         color="#17A77E", lw=2, label="Línea ideal")

# --- Línea de tendencia sobre test ---
coef_linea = np.polyfit(y_test, y_test_pred, 1)  # pendiente, intercepto
x_tendencia = np.linspace(min_val, max_val, 100)
y_tendencia = coef_linea[0] * x_tendencia + coef_linea[1]
plt.plot(x_tendencia, y_tendencia, '-', color='black',
         lw=2, label="Línea de tendencia (test)")

# --- Etiquetas ---
plt.xlabel("Turbidez real (NTU)")
plt.ylabel("Turbidez estimada (NTU)")

plt.title(
    f"Modelo XGBoost (Optimizado con Optuna)\n"
    f"R² train: {r2_train:.3f}, RMSE: {rmse_train:.3f} | "
    f"R² test: {r2_test:.3f}, RMSE: {rmse_test:.3f}"
)

plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()
C:\Users\VICTOR\AppData\Local\Programs\Python\Python313\Lib\site-packages\tqdm\auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm
  0%|          | 0/50 [00:00<?, ?it/s]Best trial: 0. Best value: 0.897072:   0%|          | 0/50 [00:00<?, ?it/s]Best trial: 0. Best value: 0.897072:   2%|▏         | 1/50 [00:00<00:17,  2.76it/s]Best trial: 0. Best value: 0.897072:   2%|▏         | 1/50 [00:00<00:17,  2.76it/s]Best trial: 0. Best value: 0.897072:   4%|▍         | 2/50 [00:00<00:14,  3.25it/s]Best trial: 2. Best value: 0.922563:   4%|▍         | 2/50 [00:01<00:14,  3.25it/s]Best trial: 2. Best value: 0.922563:   6%|▌         | 3/50 [00:01<00:16,  2.91it/s]Best trial: 2. Best value: 0.922563:   6%|▌         | 3/50 [00:01<00:16,  2.91it/s]Best trial: 2. Best value: 0.922563:   8%|▊         | 4/50 [00:01<00:14,  3.18it/s]Best trial: 2. Best value: 0.922563:   8%|▊         | 4/50 [00:01<00:14,  3.18it/s]Best trial: 2. Best value: 0.922563:  10%|█         | 5/50 [00:01<00:17,  2.53it/s]Best trial: 2. Best value: 0.922563:  10%|█         | 5/50 [00:02<00:17,  2.53it/s]Best trial: 2. Best value: 0.922563:  12%|█▏        | 6/50 [00:02<00:20,  2.12it/s]Best trial: 2. Best value: 0.922563:  12%|█▏        | 6/50 [00:03<00:20,  2.12it/s]Best trial: 2. Best value: 0.922563:  14%|█▍        | 7/50 [00:03<00:21,  1.96it/s]Best trial: 2. Best value: 0.922563:  14%|█▍        | 7/50 [00:03<00:21,  1.96it/s]Best trial: 2. Best value: 0.922563:  16%|█▌        | 8/50 [00:03<00:20,  2.02it/s]Best trial: 2. Best value: 0.922563:  16%|█▌        | 8/50 [00:03<00:20,  2.02it/s]Best trial: 2. Best value: 0.922563:  18%|█▊        | 9/50 [00:03<00:17,  2.36it/s]Best trial: 2. Best value: 0.922563:  18%|█▊        | 9/50 [00:04<00:17,  2.36it/s]Best trial: 2. Best value: 0.922563:  20%|██        | 10/50 [00:04<00:19,  2.08it/s]Best trial: 2. Best value: 0.922563:  20%|██        | 10/50 [00:05<00:19,  2.08it/s]Best trial: 2. Best value: 0.922563:  22%|██▏       | 11/50 [00:05<00:20,  1.88it/s]Best trial: 2. Best value: 0.922563:  22%|██▏       | 11/50 [00:05<00:20,  1.88it/s]Best trial: 2. Best value: 0.922563:  24%|██▍       | 12/50 [00:05<00:21,  1.75it/s]Best trial: 2. Best value: 0.922563:  24%|██▍       | 12/50 [00:06<00:21,  1.75it/s]Best trial: 2. Best value: 0.922563:  26%|██▌       | 13/50 [00:06<00:22,  1.64it/s]Best trial: 2. Best value: 0.922563:  26%|██▌       | 13/50 [00:07<00:22,  1.64it/s]Best trial: 2. Best value: 0.922563:  28%|██▊       | 14/50 [00:07<00:22,  1.59it/s]Best trial: 2. Best value: 0.922563:  28%|██▊       | 14/50 [00:07<00:22,  1.59it/s]Best trial: 2. Best value: 0.922563:  30%|███       | 15/50 [00:07<00:23,  1.48it/s]Best trial: 2. Best value: 0.922563:  30%|███       | 15/50 [00:08<00:23,  1.48it/s]Best trial: 2. Best value: 0.922563:  32%|███▏      | 16/50 [00:08<00:21,  1.57it/s]Best trial: 16. Best value: 0.923312:  32%|███▏      | 16/50 [00:08<00:21,  1.57it/s]Best trial: 16. Best value: 0.923312:  34%|███▍      | 17/50 [00:08<00:19,  1.67it/s]Best trial: 17. Best value: 0.933624:  34%|███▍      | 17/50 [00:09<00:19,  1.67it/s]Best trial: 17. Best value: 0.933624:  36%|███▌      | 18/50 [00:09<00:17,  1.82it/s]Best trial: 17. Best value: 0.933624:  36%|███▌      | 18/50 [00:09<00:17,  1.82it/s]Best trial: 17. Best value: 0.933624:  38%|███▊      | 19/50 [00:09<00:15,  2.03it/s]Best trial: 17. Best value: 0.933624:  38%|███▊      | 19/50 [00:10<00:15,  2.03it/s]Best trial: 17. Best value: 0.933624:  40%|████      | 20/50 [00:10<00:15,  1.96it/s]Best trial: 17. Best value: 0.933624:  40%|████      | 20/50 [00:10<00:15,  1.96it/s]Best trial: 17. Best value: 0.933624:  42%|████▏     | 21/50 [00:10<00:13,  2.11it/s]Best trial: 17. Best value: 0.933624:  42%|████▏     | 21/50 [00:10<00:13,  2.11it/s]Best trial: 17. Best value: 0.933624:  44%|████▍     | 22/50 [00:10<00:11,  2.39it/s]Best trial: 17. Best value: 0.933624:  44%|████▍     | 22/50 [00:11<00:11,  2.39it/s]Best trial: 17. Best value: 0.933624:  46%|████▌     | 23/50 [00:11<00:10,  2.46it/s]Best trial: 17. Best value: 0.933624:  46%|████▌     | 23/50 [00:11<00:10,  2.46it/s]Best trial: 17. Best value: 0.933624:  48%|████▊     | 24/50 [00:11<00:11,  2.34it/s]Best trial: 17. Best value: 0.933624:  48%|████▊     | 24/50 [00:12<00:11,  2.34it/s]Best trial: 17. Best value: 0.933624:  50%|█████     | 25/50 [00:12<00:10,  2.36it/s]Best trial: 17. Best value: 0.933624:  50%|█████     | 25/50 [00:12<00:10,  2.36it/s]Best trial: 17. Best value: 0.933624:  52%|█████▏    | 26/50 [00:12<00:11,  2.09it/s]Best trial: 17. Best value: 0.933624:  52%|█████▏    | 26/50 [00:13<00:11,  2.09it/s]Best trial: 17. Best value: 0.933624:  54%|█████▍    | 27/50 [00:13<00:11,  1.92it/s]Best trial: 17. Best value: 0.933624:  54%|█████▍    | 27/50 [00:13<00:11,  1.92it/s]Best trial: 17. Best value: 0.933624:  56%|█████▌    | 28/50 [00:13<00:11,  1.86it/s]Best trial: 17. Best value: 0.933624:  56%|█████▌    | 28/50 [00:14<00:11,  1.86it/s]Best trial: 17. Best value: 0.933624:  58%|█████▊    | 29/50 [00:14<00:11,  1.91it/s]Best trial: 17. Best value: 0.933624:  58%|█████▊    | 29/50 [00:14<00:11,  1.91it/s]Best trial: 17. Best value: 0.933624:  60%|██████    | 30/50 [00:14<00:10,  2.00it/s]Best trial: 17. Best value: 0.933624:  60%|██████    | 30/50 [00:15<00:10,  2.00it/s]Best trial: 17. Best value: 0.933624:  62%|██████▏   | 31/50 [00:15<00:10,  1.85it/s]Best trial: 17. Best value: 0.933624:  62%|██████▏   | 31/50 [00:16<00:10,  1.85it/s]Best trial: 17. Best value: 0.933624:  64%|██████▍   | 32/50 [00:16<00:09,  1.92it/s]Best trial: 17. Best value: 0.933624:  64%|██████▍   | 32/50 [00:16<00:09,  1.92it/s]Best trial: 17. Best value: 0.933624:  66%|██████▌   | 33/50 [00:16<00:09,  1.87it/s]Best trial: 17. Best value: 0.933624:  66%|██████▌   | 33/50 [00:17<00:09,  1.87it/s]Best trial: 17. Best value: 0.933624:  68%|██████▊   | 34/50 [00:17<00:08,  1.90it/s]Best trial: 17. Best value: 0.933624:  68%|██████▊   | 34/50 [00:17<00:08,  1.90it/s]Best trial: 17. Best value: 0.933624:  70%|███████   | 35/50 [00:17<00:07,  1.93it/s]Best trial: 17. Best value: 0.933624:  70%|███████   | 35/50 [00:18<00:07,  1.93it/s]Best trial: 17. Best value: 0.933624:  72%|███████▏  | 36/50 [00:18<00:07,  1.81it/s]Best trial: 17. Best value: 0.933624:  72%|███████▏  | 36/50 [00:18<00:07,  1.81it/s]Best trial: 17. Best value: 0.933624:  74%|███████▍  | 37/50 [00:18<00:06,  1.89it/s]Best trial: 17. Best value: 0.933624:  74%|███████▍  | 37/50 [00:18<00:06,  1.89it/s]Best trial: 17. Best value: 0.933624:  76%|███████▌  | 38/50 [00:18<00:05,  2.22it/s]Best trial: 17. Best value: 0.933624:  76%|███████▌  | 38/50 [00:19<00:05,  2.22it/s]Best trial: 17. Best value: 0.933624:  78%|███████▊  | 39/50 [00:19<00:05,  2.16it/s]Best trial: 17. Best value: 0.933624:  78%|███████▊  | 39/50 [00:19<00:05,  2.16it/s]Best trial: 17. Best value: 0.933624:  80%|████████  | 40/50 [00:19<00:04,  2.35it/s]Best trial: 17. Best value: 0.933624:  80%|████████  | 40/50 [00:20<00:04,  2.35it/s]Best trial: 17. Best value: 0.933624:  82%|████████▏ | 41/50 [00:20<00:03,  2.38it/s]Best trial: 17. Best value: 0.933624:  82%|████████▏ | 41/50 [00:20<00:03,  2.38it/s]Best trial: 17. Best value: 0.933624:  84%|████████▍ | 42/50 [00:20<00:03,  2.28it/s]Best trial: 17. Best value: 0.933624:  84%|████████▍ | 42/50 [00:21<00:03,  2.28it/s]Best trial: 17. Best value: 0.933624:  86%|████████▌ | 43/50 [00:21<00:03,  1.79it/s]Best trial: 17. Best value: 0.933624:  86%|████████▌ | 43/50 [00:21<00:03,  1.79it/s]Best trial: 17. Best value: 0.933624:  88%|████████▊ | 44/50 [00:21<00:03,  1.93it/s]Best trial: 17. Best value: 0.933624:  88%|████████▊ | 44/50 [00:22<00:03,  1.93it/s]Best trial: 17. Best value: 0.933624:  90%|█████████ | 45/50 [00:22<00:02,  1.97it/s]Best trial: 17. Best value: 0.933624:  90%|█████████ | 45/50 [00:23<00:02,  1.97it/s]Best trial: 17. Best value: 0.933624:  92%|█████████▏| 46/50 [00:23<00:02,  1.89it/s]Best trial: 17. Best value: 0.933624:  92%|█████████▏| 46/50 [00:23<00:02,  1.89it/s]Best trial: 17. Best value: 0.933624:  94%|█████████▍| 47/50 [00:23<00:01,  1.91it/s]Best trial: 17. Best value: 0.933624:  94%|█████████▍| 47/50 [00:24<00:01,  1.91it/s]Best trial: 17. Best value: 0.933624:  96%|█████████▌| 48/50 [00:24<00:01,  1.93it/s]Best trial: 17. Best value: 0.933624:  96%|█████████▌| 48/50 [00:24<00:01,  1.93it/s]Best trial: 17. Best value: 0.933624:  98%|█████████▊| 49/50 [00:24<00:00,  1.86it/s]Best trial: 17. Best value: 0.933624:  98%|█████████▊| 49/50 [00:25<00:00,  1.86it/s]Best trial: 17. Best value: 0.933624: 100%|██████████| 50/50 [00:25<00:00,  1.84it/s]Best trial: 17. Best value: 0.933624: 100%|██████████| 50/50 [00:25<00:00,  1.98it/s]

=== RESULTADOS OPTUNA ===
Mejores parámetros encontrados:
max_depth learning_rate subsample colsample_bytree n_estimators R2_promedio_CV
0 6 0.071912 0.672075 0.942038 318 0.933624 
=== Resultados Entrenamiento ===
R² train: 1.000
RMSE train: 0.082

=== Resultados Testeo ===
R² test: 0.942
RMSE test: 9.568

10.6 Relevancia de Bandas espectrales en el modelo XGBOOST

La importancia de las bandas indica cuánto contribuye cada banda espectral (o combinación de ellas) al rendimiento del modelo. Esto permite saber cuáles longitudes de onda son más informativas para estimar la variable objetivo.

Interpretación:

Cuanto mayor sea la importancia, más influye esa banda en las divisiones del modelo.

No necesariamente significa correlación directa (puede ser no lineal).

Si una banda tiene muy baja importancia (≈0), puede eliminarse en futuras versiones del modelo para simplificarlo.

Cuando entrenás un modelo como XGBoost, el algoritmo aprende qué variables (bandas espectrales) son más útiles para predecir la turbidez.

• La importancia de las bandas indica cuánto aportó cada banda al modelo para mejorar sus predicciones.

En XGBoost, el atributo model.feature_importances_ mide la importancia promedio de cada variable dentro de todos los árboles del modelo.

XGBoost da la importancia relativa (no absoluta ni porcentual) basada en “weight”. Los valores se normalizan para que la suma de todos sea 1.0.

Tipo de importancia (importance_type) Qué mide exactamente Interpretación práctica
"weight" (por defecto) Número de veces que la variable fue usada para dividir nodos Bandas usadas más veces para tomar decisiones
"gain" Promedio de la mejora en la predicción al usar esa banda Cuánto “valor” aporta cada vez que se usa
"cover" Promedio de la cantidad de datos afectados por esas divisiones Qué tan general es su influencia
"total_gain" o "total_cover" Suma total de las métricas anteriores Importancia acumulada

Bandas de Sentinel-2

Bandas de Sentinel-2

Importancia de cada banda dentro del modelo XGBoost

Código 
# === IMPORTANCIA DE BANDAS (Feature Importances) ===
import pandas as pd
import matplotlib.pyplot as plt

# Obtener importancias del modelo entrenado
importancias = model.feature_importances_

# Crear DataFrame ordenado
importancia_df = pd.DataFrame({
    'Banda': X.columns,
    'Importancia': importancias
}).sort_values(by='Importancia', ascending=False)

display(importancia_df)

# Gráfico de barras ordenado
plt.figure(figsize=(8, 4))
plt.bar(importancia_df['Banda'], importancia_df['Importancia'], color='#2E86AB')
plt.title('Importancia de las bandas en el modelo XGBoost')
plt.xlabel('Banda')
plt.ylabel('Importancia relativa')
plt.xticks(rotation=45)
plt.grid(axis='y', linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()
Banda Importancia
4 B05 0.719612
5 B06 0.157507
0 B01 0.052247
7 B08 0.027199
1 B02 0.017780
10 B8A 0.009735
8 B11 0.005484
6 B07 0.002954
3 B04 0.002905
2 B03 0.002719
9 B12 0.001858

B05 explica la mayor parte del poder predictivo total del modelo.

Las bandas con valores cercanos a cero aportan muy poca información o son redundantes con otras.

Las bandas Red Edge (B05–B06-B07) se ubican justo en la transición entre el rojo y el infrarrojo cercano, donde la reflectancia del agua cambia mucho con la concentración de sedimentos. → Son excelentes indicadores de turbidez o material particulado suspendido.

Las bandas azules (B01–B02) todavía aportan algo, porque la turbidez afecta la dispersión del azul. → En aguas claras, el azul predomina; en aguas turbias, se atenúa.

Las bandas SWIR (B11, B12) y el rojo visible (B04) no aportan casi nada, ya que esas longitudes de onda no penetran el agua y son más útiles para suelo o vegetación.

10.7 Mapas XGBOOST con hiperparámetros optimizados

Código 
import rasterio as rs
import numpy as np
import matplotlib.pyplot as plt
import glob
import os
import pandas as pd
from xgboost import XGBRegressor
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import r2_score, mean_squared_error
import optuna
from skimage.filters import threshold_otsu

# --- CONFIGURACIÓN ---
carpeta = r"D:\GIT\Turbidez\recorte_acolite"
salida = os.path.join(carpeta, "turbidez_out_XGBOOST_optuna")
os.makedirs(salida, exist_ok=True)

# --- ENTRENAMIENTO DEL MODELO OPTIMIZADO ---
df = pd.read_csv(r"D:\GIT\Turbidez\Datos creados\Datos_turb_banda\C1_turb_banda.csv")

# Variables predictoras y objetivo
X = df[["B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08", "B11", "B12", "B8A"]]
y = df["turbidez"]

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=42)
optuna.logging.set_verbosity(optuna.logging.WARNING)

def objective(trial):
    params = {
        'max_depth': trial.suggest_int('max_depth', 3, 10),
        'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.3),
        'subsample': trial.suggest_float('subsample', 0.6, 1.0),
        'colsample_bytree': trial.suggest_float('colsample_bytree', 0.6, 1.0),
        'n_estimators': trial.suggest_int('n_estimators', 200, 800),
        'random_state': 42
    }
    model = XGBRegressor(**params)
    score = cross_val_score(model, X_train, y_train, scoring='r2', cv=5).mean()
    return score

study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=50, show_progress_bar=False)

best_params = study.best_params
print("=== MEJORES PARÁMETROS ===")
print(best_params)
print(f"R² promedio CV: {study.best_value:.3f}")

# Entrenar modelo final
model = XGBRegressor(**best_params)
model.fit(X_train, y_train)

# Evaluar desempeño
y_test_pred = model.predict(X_test)
print(f"R² test: {r2_score(y_test, y_test_pred):.3f}")
print(f"RMSE test: {np.sqrt(mean_squared_error(y_test, y_test_pred)):.3f}")

# --- PROCESAR RASTERS ---
feature_names = ["B01", "B02", "B03", "B04", "B05",
                 "B06", "B07", "B08", "B11", "B12", "B8A"]

tifs = sorted(glob.glob(os.path.join(carpeta, "*.tif")))
print(f"\nSe encontraron {len(tifs)} archivos .tif")

for ruta in tifs:
    nombre = os.path.basename(ruta)
    fecha = os.path.splitext(nombre)[0]
    print(f"\nProcesando: {nombre}")

    with rs.open(ruta) as raster:
        bandas = [raster.read(i).astype(np.float32) for i in range(1, 12)]
        perfil = raster.profile

    # --- NDWI y máscara de agua ---
    B03, B11 = bandas[2], bandas[9]
    ndwi = (B03 - B11) / (B03 + B11 + 1e-6)
    umbral = np.nanmean(ndwi)
    mascara = (ndwi > umbral).astype(np.uint8)

    # --- Crear matriz de predicción ---
    filas, cols = B03.shape
    X_img = np.stack(bandas, axis=-1)
    X_img_2d = X_img.reshape(-1, X_img.shape[-1])
    X_img_df = pd.DataFrame(X_img_2d, columns=feature_names)

    # --- Predecir turbidez ---
    turbidez_pred = model.predict(X_img_df)
    turbidez_pred = turbidez_pred.reshape(filas, cols)

    # Aplicar máscara
    turbidez_agua = np.where(mascara == 1, turbidez_pred, np.nan)

    # --- Exportar GeoTIFF ---
    salida_tif = os.path.join(salida, f"turbidez_{fecha}.tif")
    perfil.update(dtype=rs.float32, count=1, compress='lzw', nodata=np.nan)
    with rs.open(salida_tif, 'w', **perfil) as dst:
        dst.write(turbidez_agua.astype(np.float32), 1)

    # --- Visualización ---
    p5, p95 = np.nanpercentile(turbidez_agua, [5, 95])
    fig, axes = plt.subplots(1, 2, figsize=(8, 4.5))

    im1 = axes[0].imshow(mascara, cmap="Greys")
    axes[0].set_title(f"Máscara de agua - {fecha}")
    axes[0].axis("off")
    fig.colorbar(im1, ax=axes[0], fraction=0.046, pad=0.04, label="Agua / No Agua")

    im2 = axes[1].imshow(turbidez_agua, cmap="rainbow", vmin=p5, vmax=p95)
    axes[1].set_title(f"Turbidez (XGBoost Optuna) - {fecha}")
    axes[1].axis("off")
    fig.colorbar(im2, ax=axes[1], fraction=0.046, pad=0.04, label="Turbidez (NTU)")

    plt.tight_layout()
    plt.show()
=== MEJORES PARÁMETROS ===
{'max_depth': 3, 'learning_rate': 0.022172460819822562, 'subsample': 0.7605112208446342, 'colsample_bytree': 0.7651599773171959, 'n_estimators': 372}
R² promedio CV: 0.935
R² test: 0.902
RMSE test: 12.482

Se encontraron 7 archivos .tif

Procesando: 2023-05-11.tif


Procesando: 2023-12-12.tif


Procesando: 2024-05-20.tif


Procesando: 2024-05-30.tif


Procesando: 2024-12-16.tif


Procesando: 2025-06-09.tif


Procesando: 2025-09-12.tif

Referencias

[1]
J. cc Delegido et al., «Turbidity and secchi disc depth with sentinel-2 in different trophic status reservoirs at the comunidad valenciana», Revista de Teledeteccion, vol. 2019, pp. 15-24, 2019, doi: 10.4995/raet.2019.12603.
[2]
Y. Ma et al., «Remote sensing of turbidity for lakes in Northeast China using sentinel-2 images with machine learning algorithms», IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 14, pp. 9132-9146, 2021, doi: 10.1109/JSTARS.2021.3109292.
[3]
A. K. M. A. Hossain, C. Mathias, y R. Blanton, «Remote sensing of turbidity in the tennessee river using landsat 8 satellite», Remote Sensing, vol. 13, sep. 2021, doi: 10.3390/rs13183785.
[4]
M. Najafzadeh y S. Basirian, «Evaluation of River Water Quality Index Using Remote Sensing and Artificial Intelligence Models», Remote Sensing, vol. 15, may 2023, doi: 10.3390/rs15092359.
[5]
M. Allam, M. Y. A. Khan, y Q. Meng, «Retrieval of turbidity on a spatio-temporal scale using landsat 8 SR: A case study of the Ramganga River in the Ganges Basin, India», Applied Sciences (Switzerland), vol. 10, jun. 2020, doi: 10.3390/app10113702.
[6]
Y. Luo, D. Doxaran, y Q. Vanhellemont, «Retrieval and validation of water turbidity at metre-scale using Pleiades satellite data: A case study in the Gironde Estuary», Remote Sensing, vol. 12, mar. 2020, doi: 10.3390/rs12060946.
[7]
M. Shen, S. Wang, Y. Li, M. Tang, y Y. Ma, «Pattern of turbidity change in the middle reaches of the yarlung zangbo river, southern tibetan plateau, from 2007 to 2017», Remote Sensing, vol. 13, pp. 1-25, ene. 2021, doi: 10.3390/rs13020182.
[8]
Y. O. Ouma, K. Noor, y K. Herbert, «Modelling Reservoir Chlorophyll- a, TSS, and Turbidity Using Sentinel-2A MSI and Landsat-8 OLI Satellite Sensors with Empirical Multivariate Regression», Journal of Sensors, vol. 2020, 2020, doi: 10.1155/2020/8858408.
[9]
Z. Yang et al., «Combined Retrievals of Turbidity from Sentinel-2A/B and Landsat-8/9 in the Taihu Lake through Machine Learning», Remote Sensing, vol. 15, sep. 2023, doi: 10.3390/rs15174333.
[10]
Y. Zhao et al., «Retrieval of Water Quality Parameters Based on Near-Surface Remote Sensing and Machine Learning Algorithm», Remote Sensing, vol. 14, nov. 2022, doi: 10.3390/rs14215305.
[11]
S. Chatterjee y A. S. Hadi, Regression analysis by example. John Wiley & Sons, 2015.

Notas

  1. 0,83 - 112,26 NTU.↩︎

  2. 2,3 - 107,02 NTU.↩︎

  3. IOA = index of agreement
    SI = scatter index.↩︎

  4. 20,6 - 112 NTU
    2,3 - 15,4 NTU.↩︎

  5. MAPE = Mean Absolute Percentage Error
    0 - 1300 NTU
    0 - 80 NTU.↩︎